From 6cdebb3f6d5d65f699d8fe00ef793806f48ea86a Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 9 May 2025 17:50:11 +0100 Subject: [PATCH 1/5] added Jack's dPIE --- autogalaxy/profiles/mass/total/__init__.py | 1 + .../mass/total/dual_pseudo_isothermal.py | 221 ++++++++++++++++++ 2 files changed, 222 insertions(+) create mode 100644 autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py diff --git a/autogalaxy/profiles/mass/total/__init__.py b/autogalaxy/profiles/mass/total/__init__.py index 9d7b096ad..430ec2674 100644 --- a/autogalaxy/profiles/mass/total/__init__.py +++ b/autogalaxy/profiles/mass/total/__init__.py @@ -1,3 +1,4 @@ +from .dual_pseudo_isothermal import dPIE, dPIESph 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.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py new file mode 100644 index 000000000..ca9cd7261 --- /dev/null +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py @@ -0,0 +1,221 @@ +from typing import Tuple +import numpy as np + +import autoarray as aa +from autogalaxy.profiles.mass.abstract.abstract import MassProfile + + +class dPIESph(MassProfile): + ''' + The dual Pseudo-Isothermal Elliptical mass distribution introduced in + Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + + This version is without ellipticity, so perhaps the "E" is a misnomer. + + Corresponds to a projected density profile that looks like: + + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + + (c.f. Eliasdottir '07 eqn. A3) + + In this parameterization, ra and rs are the scale radii above in angular + units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. + ''' + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + ra: float = 0.1, + rs: float = 2.0, + kappa_scale: float = 0.1, + ): + # Ensure rs > ra (things will probably break otherwise) + if ra > rs: + ra, rs = rs, ra + super(MassProfile, self).__init__(centre, (0.0, 0.0)) + self.ra = ra + self.rs = rs + self.kappa_scale = kappa_scale + + def _deflection_angle(self, radii): + ''' + For a circularly symmetric dPIE profile, computes the magnitude of the deflection at each radius. + ''' + r_ra = radii / self.ra + r_rs = radii / self.rs + # c.f. Eliasdottir '07 eq. A20 + f = ( + r_ra / (1 + np.sqrt(1 + r_ra * r_ra)) + - r_rs / (1 + np.sqrt(1 + r_rs * r_rs)) + ) + + ra, rs = self.ra, self.rs + # c.f. Eliasdottir '07 eq. A19 + # magnitude of deflection + alpha = 2 * self.kappa_scale * ra * rs / (rs - ra) * f + return alpha + + def _convergence(self, radii): + radsq = radii * radii + a, s = self.ra, self.rs + # c.f. Eliasdottir '07 eqn (A3) + return ( + self.kappa_scale * (a * s) / (s - a) * + (1/np.sqrt(a**2 + radsq) - 1/np.sqrt(s**2 + radsq)) + ) + + def _potential(self, radii): + raise NotImplementedError + + @aa.grid_dec.to_vector_yx + def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + ys, xs = grid.T + (ycen, xcen) = self.centre + xoff, yoff = xs - xcen, ys - ycen + radii = np.sqrt(xoff**2 + yoff**2) + + alpha = self._deflection_angle(radii) + + # now we decompose the deflection into y/x components + defl_y = alpha * yoff / radii + defl_x = alpha * xoff / radii + return aa.Grid2DIrregular.from_yx_1d( + defl_y, defl_x + ) + + @aa.grid_dec.to_array + @aa.grid_dec.transform + @aa.grid_dec.relocate_to_radial_minimum + def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + # already transformed to center on profile centre so this works + radsq = (grid[:, 0]**2 + grid[:, 1]**2) + return self._convergence(np.sqrt(radsq)) + + @aa.grid_dec.to_array + @aa.grid_dec.relocate_to_radial_minimum + def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + # already transformed to center on profile centre so this works + radsq = (grid[:, 0]**2 + grid[:, 1]**2) + return self._potential(np.sqrt(radsq)) + + +class dPIE(dPIESph): + ''' + The dual Pseudo-Isothermal Elliptical mass distribution introduced in + Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + + Corresponds to a projected density profile that looks like: + + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + + (c.f. Eliasdottir '07 eqn. A3) + + In this parameterization, ra and rs are the scale radii above in angular + units (arcsec). The parameter kappa_scale is \\Sigma_0 / \\Sigma_crit. + + WARNING: This uses the "pseud-elliptical" approximation, where the ellipticity + is applied to the *potential* rather than the *mass* to ease calculation. + Use at your own risk! (And TODO Jack: fix this!) + This approximation is used by the lenstronomy PJAFFE profile (which is the + same functional form), but not by the lenstool PIEMD (also synonymous with this), + which correctly solved the differential equations for the mass-based ellipticity. + ''' + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + ell_comps: Tuple[float, float] = (0.0, 0.0), + ra: float = 0.1, + rs: float = 2.0, + kappa_scale: float = 0.1, + ): + super(MassProfile, self).__init__(centre, ell_comps) + if ra > rs: + ra, rs = rs, ra + self.ra = ra + self.rs = rs + self.kappa_scale = kappa_scale + + def _align_to_major_axis(self, ys, xs): + ''' + Aligns coordinates to the major axis of this halo. Returns (y', x'), + where x' is along the major axis and y' is along the minor axis. + + Does NOT translate, only rotates. + ''' + costheta, sintheta = self._cos_and_sin_to_x_axis() + _xs = (costheta * xs + sintheta * ys) + _ys = (-sintheta * xs + costheta * ys) + return _ys, _xs + + def _align_from_major_axis(self, _ys, _xs): + ''' + Given _ys and _xs as offsets along the minor and major axes, + respectively, this transforms them back to the regular coordinate + system. + + Does NOT translate, only rotates. + ''' + costheta, sintheta = self._cos_and_sin_to_x_axis() + xs = (costheta * _xs + -sintheta * _ys) + ys = (sintheta * _xs + costheta * _ys) + return ys, xs + + def _ellip(self): + ellip = np.sqrt(self.ell_comps[0]**2 + self.ell_comps[1]**2) + MAX_ELLIP = 0.99999 + return min(ellip, MAX_ELLIP) + + @aa.grid_dec.to_vector_yx + def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + ys, xs = grid.T + (ycen, xcen) = self.centre + xoff, yoff = xs - xcen, ys - ycen + _ys, _xs = self._align_to_major_axis(yoff, xoff) + + ellip = self._ellip() + _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) + + # Compute the deflection magnitude of a *non-elliptical* profile + alpha_circ = self._deflection_angle(_radii) + + # This is in axes aligned to the major/minor axis + _defl_xs = alpha_circ * np.sqrt(1 - ellip) * (_xs / _radii) + _defl_ys = alpha_circ * np.sqrt(1 + ellip) * (_ys / _radii) + + # And here we convert back to the real axes + defl_ys, defl_xs = self._align_from_major_axis(_defl_ys, _defl_xs) + return aa.Grid2DIrregular.from_yx_1d( + defl_ys, defl_xs + ) + + @aa.grid_dec.to_array + def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + ys, xs = grid.T + (ycen, xcen) = self.centre + xoff, yoff = xs - xcen, ys - ycen + _ys, _xs = self._align_to_major_axis(yoff, xoff) + + ellip = self._ellip() + _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) + + # Compute the convergence and deflection of a *circular* profile + kappa_circ = self._convergence(_radii) + alpha_circ = self._deflection_angle(_radii) + + asymm_term = (ellip * (1 - ellip) * _xs**2 - ellip * (1 + ellip) * _ys**2) / _radii**2 + + # convergence = 1/2 \nabla \alpha = 1/2 \nabla^2 potential + # The "asymm_term" is asymmetric on x and y, so averages out to + # zero over all space + return kappa_circ * (1 - asymm_term) + (alpha_circ / _radii) * asymm_term + + @aa.grid_dec.to_array + def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + ys, xs = grid.T + (ycen, xcen) = self.centre + xoff, yoff = xs - xcen, ys - ycen + _ys, _xs = self._align_to_major_axis(yoff, xoff) + ellip = self._ellip() + _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) + return super(dPIESph, self)._potential(_radii) \ No newline at end of file From 41e7dd7fcfa936c1e2305401e50187874c7e2493 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 9 May 2025 19:24:15 +0100 Subject: [PATCH 2/5] added Jack's dPIE --- .../mass/total/test_dual_pseudo_isothermal.py | 175 ++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py new file mode 100644 index 000000000..28ef2b039 --- /dev/null +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py @@ -0,0 +1,175 @@ +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(): + mp = ag.mp.IsothermalSph(centre=(-0.7, 0.5), einstein_radius=1.3) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(1.21510, 1e-4) + assert deflections[0, 1] == pytest.approx(-0.46208, 1e-4) + + mp = ag.mp.IsothermalSph(centre=(-0.1, 0.1), einstein_radius=5.0) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(4.88588, 1e-4) + assert deflections[0, 1] == pytest.approx(1.06214, 1e-4) + + mp = ag.mp.Isothermal(centre=(0, 0), ell_comps=(0.0, 0.333333), einstein_radius=1.0) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(0.79421, 1e-3) + assert deflections[0, 1] == pytest.approx(0.50734, 1e-3) + + mp = ag.mp.Isothermal(centre=(0, 0), ell_comps=(0.0, 0.333333), einstein_radius=1.0) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(0.79421, 1e-3) + assert deflections[0, 1] == pytest.approx(0.50734, 1e-3) + + elliptical = ag.mp.Isothermal( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 + ) + spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) + + assert elliptical.deflections_yx_2d_from(grid=grid) == pytest.approx( + spherical.deflections_yx_2d_from(grid=grid), 1e-4 + ) + + +def test__convergence_2d_from(): + # eta = 1.0 + # kappa = 0.5 * 1.0 ** 1.0 + + mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=2.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.5 * 2.0, 1e-3) + + mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=1.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.5, 1e-3) + + mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=2.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.5 * 2.0, 1e-3) + + mp = ag.mp.Isothermal( + centre=(0.0, 0.0), ell_comps=(0.0, 0.333333), einstein_radius=1.0 + ) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.66666, 1e-3) + + elliptical = ag.mp.Isothermal( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 + ) + spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) + + assert elliptical.convergence_2d_from(grid=grid) == pytest.approx( + spherical.convergence_2d_from(grid=grid), 1e-4 + ) + + +def test__potential_2d_from(): + mp = ag.mp.IsothermalSph(centre=(-0.7, 0.5), einstein_radius=1.3) + + potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) + + assert potential == pytest.approx(1.23435, 1e-3) + + mp = ag.mp.Isothermal( + centre=(-0.7, 0.5), + ell_comps=(0.152828, -0.088235), + einstein_radius=1.3, + ) + + potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) + + assert potential == pytest.approx(1.19268, 1e-3) + + elliptical = ag.mp.Isothermal( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 + ) + spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) + + assert elliptical.potential_2d_from(grid=grid) == pytest.approx( + spherical.potential_2d_from(grid=grid), 1e-4 + ) + + +def test__shear_yx_2d_from(): + mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=2.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert shear[0, 0] == pytest.approx(0.0, 1e-4) + assert shear[0, 1] == pytest.approx(-convergence, 1e-4) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[2.0, 1.0]])) + shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[2.0, 1.0]])) + + assert shear[0, 0] == pytest.approx(-(4.0 / 5.0) * convergence, 1e-4) + assert shear[0, 1] == pytest.approx((3.0 / 5.0) * convergence, 1e-4) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[3.0, 5.0]])) + shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[3.0, 5.0]])) + + assert shear[0, 0] == pytest.approx(-(30.0 / 34.0) * convergence, 1e-4) + assert shear[0, 1] == pytest.approx(-(16.0 / 34.0) * convergence, 1e-4) + + mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=2.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert shear[0, 0] == pytest.approx(0.0, 1e-4) + assert shear[0, 1] == pytest.approx(-convergence, 1e-4) + + mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.3, 0.4), einstein_radius=2.0) + + shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert shear[0, 0] == pytest.approx(0.0, 1e-4) + assert shear[0, 1] == pytest.approx(-1.11803398874, 1e-4) + + +def test__compare_to_cored_power_law(): + isothermal = ag.mp.Isothermal( + centre=(0.0, 0.0), ell_comps=(0.333333, 0.0), einstein_radius=1.0 + ) + cored_power_law = ag.mp.PowerLawCore( + centre=(0.0, 0.0), + ell_comps=(0.333333, 0.0), + einstein_radius=1.0, + core_radius=0.0, + ) + + assert isothermal.potential_2d_from(grid=grid) == pytest.approx( + cored_power_law.potential_2d_from(grid=grid), 1e-3 + ) + assert isothermal.potential_2d_from(grid=grid) == pytest.approx( + cored_power_law.potential_2d_from(grid=grid), 1e-3 + ) + assert isothermal.deflections_yx_2d_from(grid=grid) == pytest.approx( + cored_power_law.deflections_yx_2d_from(grid=grid), 1e-3 + ) + assert isothermal.deflections_yx_2d_from(grid=grid) == pytest.approx( + cored_power_law.deflections_yx_2d_from(grid=grid), 1e-3 + ) From 845fee15a0ec73a186f73d354eb6fad4549c46ba Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 9 May 2025 19:29:00 +0100 Subject: [PATCH 3/5] DPIE plus tests added --- autogalaxy/profiles/mass/__init__.py | 2 + .../mass/total/dual_pseudo_isothermal.py | 59 ++-- .../mass/total/test_dual_pseudo_isothermal.py | 260 ++++++------------ 3 files changed, 117 insertions(+), 204 deletions(-) diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index ab248dffd..5080042a4 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -1,6 +1,8 @@ from .abstract.abstract import MassProfile from .point import PointMass, SMBH, SMBHBinary from .total import ( + dPIE, + dPIESph, PowerLawCore, PowerLawCoreSph, PowerLawBroken, diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py index ca9cd7261..e4ee7c2cc 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py @@ -6,22 +6,7 @@ class dPIESph(MassProfile): - ''' - The dual Pseudo-Isothermal Elliptical mass distribution introduced in - Eliasdottir 2007: https://arxiv.org/abs/0710.5636 - - This version is without ellipticity, so perhaps the "E" is a misnomer. - Corresponds to a projected density profile that looks like: - - \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * - (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) - - (c.f. Eliasdottir '07 eqn. A3) - - In this parameterization, ra and rs are the scale radii above in angular - units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. - ''' def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -29,6 +14,34 @@ def __init__( rs: float = 2.0, kappa_scale: float = 0.1, ): + """ + The dual Pseudo-Isothermal Elliptical mass distribution introduced in + Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + + This version is without ellipticity, so perhaps the "E" is a misnomer. + + Corresponds to a projected density profile that looks like: + + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + + (c.f. Eliasdottir '07 eqn. A3) + + In this parameterization, ra and rs are the scale radii above in angular + units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + ra + A scale radius in arc-seconds. + rs + The second scale radius in arc-seconds. + kappa_scale + Scales the overall normalization of the profile, so related to the mass + """ + # Ensure rs > ra (things will probably break otherwise) if ra > rs: ra, rs = rs, ra @@ -64,9 +77,6 @@ def _convergence(self, radii): (1/np.sqrt(a**2 + radsq) - 1/np.sqrt(s**2 + radsq)) ) - def _potential(self, radii): - raise NotImplementedError - @aa.grid_dec.to_vector_yx def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): ys, xs = grid.T @@ -92,11 +102,8 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): return self._convergence(np.sqrt(radsq)) @aa.grid_dec.to_array - @aa.grid_dec.relocate_to_radial_minimum def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - # already transformed to center on profile centre so this works - radsq = (grid[:, 0]**2 + grid[:, 1]**2) - return self._potential(np.sqrt(radsq)) + return np.zeros(shape=grid.shape[0]) class dPIE(dPIESph): @@ -212,10 +219,4 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @aa.grid_dec.to_array def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - ys, xs = grid.T - (ycen, xcen) = self.centre - xoff, yoff = xs - xcen, ys - ycen - _ys, _xs = self._align_to_major_axis(yoff, xoff) - ellip = self._ellip() - _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) - return super(dPIESph, self)._potential(_radii) \ No newline at end of file + return np.zeros(shape=grid.shape[0]) \ No newline at end of file diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py index 28ef2b039..988fe2ad8 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py @@ -1,175 +1,85 @@ -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(): - mp = ag.mp.IsothermalSph(centre=(-0.7, 0.5), einstein_radius=1.3) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(1.21510, 1e-4) - assert deflections[0, 1] == pytest.approx(-0.46208, 1e-4) - - mp = ag.mp.IsothermalSph(centre=(-0.1, 0.1), einstein_radius=5.0) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(4.88588, 1e-4) - assert deflections[0, 1] == pytest.approx(1.06214, 1e-4) - - mp = ag.mp.Isothermal(centre=(0, 0), ell_comps=(0.0, 0.333333), einstein_radius=1.0) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(0.79421, 1e-3) - assert deflections[0, 1] == pytest.approx(0.50734, 1e-3) - - mp = ag.mp.Isothermal(centre=(0, 0), ell_comps=(0.0, 0.333333), einstein_radius=1.0) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(0.79421, 1e-3) - assert deflections[0, 1] == pytest.approx(0.50734, 1e-3) - - elliptical = ag.mp.Isothermal( - centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 - ) - spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) - - assert elliptical.deflections_yx_2d_from(grid=grid) == pytest.approx( - spherical.deflections_yx_2d_from(grid=grid), 1e-4 - ) - - -def test__convergence_2d_from(): - # eta = 1.0 - # kappa = 0.5 * 1.0 ** 1.0 - - mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=2.0) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert convergence == pytest.approx(0.5 * 2.0, 1e-3) - - mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=1.0) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert convergence == pytest.approx(0.5, 1e-3) - - mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=2.0) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert convergence == pytest.approx(0.5 * 2.0, 1e-3) - - mp = ag.mp.Isothermal( - centre=(0.0, 0.0), ell_comps=(0.0, 0.333333), einstein_radius=1.0 - ) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert convergence == pytest.approx(0.66666, 1e-3) - - elliptical = ag.mp.Isothermal( - centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 - ) - spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) - - assert elliptical.convergence_2d_from(grid=grid) == pytest.approx( - spherical.convergence_2d_from(grid=grid), 1e-4 - ) - - -def test__potential_2d_from(): - mp = ag.mp.IsothermalSph(centre=(-0.7, 0.5), einstein_radius=1.3) - - potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) - - assert potential == pytest.approx(1.23435, 1e-3) - - mp = ag.mp.Isothermal( - centre=(-0.7, 0.5), - ell_comps=(0.152828, -0.088235), - einstein_radius=1.3, - ) - - potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) - - assert potential == pytest.approx(1.19268, 1e-3) - - elliptical = ag.mp.Isothermal( - centre=(1.1, 1.1), ell_comps=(0.0, 0.0), einstein_radius=3.0 - ) - spherical = ag.mp.IsothermalSph(centre=(1.1, 1.1), einstein_radius=3.0) - - assert elliptical.potential_2d_from(grid=grid) == pytest.approx( - spherical.potential_2d_from(grid=grid), 1e-4 - ) - - -def test__shear_yx_2d_from(): - mp = ag.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=2.0) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert shear[0, 0] == pytest.approx(0.0, 1e-4) - assert shear[0, 1] == pytest.approx(-convergence, 1e-4) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[2.0, 1.0]])) - shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[2.0, 1.0]])) - - assert shear[0, 0] == pytest.approx(-(4.0 / 5.0) * convergence, 1e-4) - assert shear[0, 1] == pytest.approx((3.0 / 5.0) * convergence, 1e-4) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[3.0, 5.0]])) - shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[3.0, 5.0]])) - - assert shear[0, 0] == pytest.approx(-(30.0 / 34.0) * convergence, 1e-4) - assert shear[0, 1] == pytest.approx(-(16.0 / 34.0) * convergence, 1e-4) - - mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), einstein_radius=2.0) - - convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert shear[0, 0] == pytest.approx(0.0, 1e-4) - assert shear[0, 1] == pytest.approx(-convergence, 1e-4) - - mp = ag.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.3, 0.4), einstein_radius=2.0) - - shear = mp.shear_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) - - assert shear[0, 0] == pytest.approx(0.0, 1e-4) - assert shear[0, 1] == pytest.approx(-1.11803398874, 1e-4) - - -def test__compare_to_cored_power_law(): - isothermal = ag.mp.Isothermal( - centre=(0.0, 0.0), ell_comps=(0.333333, 0.0), einstein_radius=1.0 - ) - cored_power_law = ag.mp.PowerLawCore( - centre=(0.0, 0.0), - ell_comps=(0.333333, 0.0), - einstein_radius=1.0, - core_radius=0.0, - ) - - assert isothermal.potential_2d_from(grid=grid) == pytest.approx( - cored_power_law.potential_2d_from(grid=grid), 1e-3 - ) - assert isothermal.potential_2d_from(grid=grid) == pytest.approx( - cored_power_law.potential_2d_from(grid=grid), 1e-3 - ) - assert isothermal.deflections_yx_2d_from(grid=grid) == pytest.approx( - cored_power_law.deflections_yx_2d_from(grid=grid), 1e-3 - ) - assert isothermal.deflections_yx_2d_from(grid=grid) == pytest.approx( - cored_power_law.deflections_yx_2d_from(grid=grid), 1e-3 - ) +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(): + mp = ag.mp.dPIESph(centre=(-0.7, 0.5), kappa_scale=1.3, 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) + + mp = ag.mp.dPIESph(centre=(-0.1, 0.1), kappa_scale=5.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) + + mp = ag.mp.dPIE(centre=(0, 0), ell_comps=(0.0, 0.333333), kappa_scale=1.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.186341843, 1e-3) + assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) + + mp = ag.mp.dPIE(centre=(0, 0), ell_comps=(0.0, 0.333333), kappa_scale=1.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.186341843, 1e-3) + assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) + + elliptical = ag.mp.dPIE( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), kappa_scale=3.0, ra=2.0, rs=3.0 + ) + spherical = ag.mp.dPIESph(centre=(1.1, 1.1), kappa_scale=3.0, ra=2.0, rs=3.0) + + assert elliptical.deflections_yx_2d_from(grid=grid) == pytest.approx( + spherical.deflections_yx_2d_from(grid=grid), 1e-4 + ) + + +def test__convergence_2d_from(): + # eta = 1.0 + # kappa = 0.5 * 1.0 ** 1.0 + + mp = ag.mp.dPIESph(centre=(0.0, 0.0), kappa_scale=2.0, ra=2.0, rs=3.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(1.57182995, 1e-3) + + mp = ag.mp.dPIE(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), kappa_scale=1.0, ra=2.0, rs=3.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.78591498, 1e-3) + + mp = ag.mp.dPIE(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), kappa_scale=2.0, ra=2.0, rs=3.0) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(1.57182995, 1e-3) + + mp = ag.mp.dPIE( + centre=(0.0, 0.0), ell_comps=(0.0, 0.333333), kappa_scale=1.0, ra=2.0, rs=3.0 + ) + + convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) + + assert convergence == pytest.approx(0.87182837, 1e-3) + + elliptical = ag.mp.dPIE( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), kappa_scale=3.0, ra=2.0, rs=3.0 + ) + spherical = ag.mp.dPIESph(centre=(1.1, 1.1), kappa_scale=3.0, ra=2.0, rs=3.0) + + assert elliptical.convergence_2d_from(grid=grid) == pytest.approx( + spherical.convergence_2d_from(grid=grid), 1e-4 + ) + From e5e562594976f13960ab2c4b7d3215cadaa5a0b3 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 13 May 2025 20:22:53 +0100 Subject: [PATCH 4/5] refactored dPIE to use autolens grid API --- .../mass/total/dual_pseudo_isothermal.py | 235 ++++++++---------- docs/installation/conda.rst | 2 +- docs/installation/pip.rst | 2 +- 3 files changed, 109 insertions(+), 130 deletions(-) diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py index e4ee7c2cc..949d1f9b0 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py @@ -5,51 +5,51 @@ from autogalaxy.profiles.mass.abstract.abstract import MassProfile -class dPIESph(MassProfile): +class dPIE(MassProfile): + ''' + The dual Pseudo-Isothermal Elliptical mass distribution introduced in + Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + + Corresponds to a projected density profile that looks like: + + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + + (c.f. Eliasdottir '07 eqn. A3) + + In this parameterization, ra and rs are the scale radii above in angular + units (arcsec). The parameter kappa_scale is \\Sigma_0 / \\Sigma_crit. + + WARNING: This uses the "pseud-elliptical" approximation, where the ellipticity + is applied to the *potential* rather than the *mass* to ease calculation. + Use at your own risk! (And TODO Jack: fix this!) + This approximation is used by the lenstronomy PJAFFE profile (which is the + same functional form), but not by the lenstool PIEMD (also synonymous with this), + which correctly solved the differential equations for the mass-based ellipticity. + ''' def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), + ell_comps: Tuple[float, float] = (0.0, 0.0), ra: float = 0.1, rs: float = 2.0, kappa_scale: float = 0.1, ): - """ - The dual Pseudo-Isothermal Elliptical mass distribution introduced in - Eliasdottir 2007: https://arxiv.org/abs/0710.5636 - - This version is without ellipticity, so perhaps the "E" is a misnomer. - - Corresponds to a projected density profile that looks like: - - \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * - (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) - - (c.f. Eliasdottir '07 eqn. A3) - - In this parameterization, ra and rs are the scale radii above in angular - units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. - - Parameters - ---------- - centre - The (y,x) arc-second coordinates of the profile centre. - ra - A scale radius in arc-seconds. - rs - The second scale radius in arc-seconds. - kappa_scale - Scales the overall normalization of the profile, so related to the mass - """ + super().__init__(centre, ell_comps) - # Ensure rs > ra (things will probably break otherwise) if ra > rs: ra, rs = rs, ra - super(MassProfile, self).__init__(centre, (0.0, 0.0)) + self.ra = ra self.rs = rs self.kappa_scale = kappa_scale + def _ellip(self): + ellip = np.sqrt(self.ell_comps[0]**2 + self.ell_comps[1]**2) + MAX_ELLIP = 0.99999 + return min(ellip, MAX_ELLIP) + def _deflection_angle(self, radii): ''' For a circularly symmetric dPIE profile, computes the magnitude of the deflection at each radius. @@ -78,144 +78,123 @@ def _convergence(self, radii): ) @aa.grid_dec.to_vector_yx + @aa.grid_dec.transform + @aa.grid_dec.relocate_to_radial_minimum def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - ys, xs = grid.T - (ycen, xcen) = self.centre - xoff, yoff = xs - xcen, ys - ycen - radii = np.sqrt(xoff**2 + yoff**2) - alpha = self._deflection_angle(radii) + ellip = self._ellip() + grid_radii = np.sqrt(grid[:,1]**2 * (1 - ellip) + grid[:,0]**2 * (1 + ellip)) - # now we decompose the deflection into y/x components - defl_y = alpha * yoff / radii - defl_x = alpha * xoff / radii - return aa.Grid2DIrregular.from_yx_1d( - defl_y, defl_x + # Compute the deflection magnitude of a *non-elliptical* profile + alpha_circ = self._deflection_angle(grid_radii) + + # This is in axes aligned to the major/minor axis + deflection_y = alpha_circ * np.sqrt(1 + ellip) * (grid[:,0] / grid_radii) + deflection_x = alpha_circ * np.sqrt(1 - ellip) * (grid[:,1] / grid_radii) + + # And here we convert back to the real axes + return self.rotated_grid_from_reference_frame_from( + grid=np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T), + **kwargs ) - @aa.grid_dec.to_array + @aa.grid_dec.to_vector_yx @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - # already transformed to center on profile centre so this works - radsq = (grid[:, 0]**2 + grid[:, 1]**2) - return self._convergence(np.sqrt(radsq)) - @aa.grid_dec.to_array - def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - return np.zeros(shape=grid.shape[0]) + ellip = self._ellip() + grid_radii = np.sqrt(grid[:,1]**2 * (1 - ellip) + grid[:,0]**2 * (1 + ellip)) + # Compute the convergence and deflection of a *circular* profile + kappa_circ = self._convergence(grid_radii) + alpha_circ = self._deflection_angle(grid_radii) -class dPIE(dPIESph): - ''' - The dual Pseudo-Isothermal Elliptical mass distribution introduced in - Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + asymm_term = (ellip * (1 - ellip) * grid[:,1]**2 - ellip * (1 + ellip) * grid[:,0]**2) / grid_radii**2 - Corresponds to a projected density profile that looks like: + # convergence = 1/2 \nabla \alpha = 1/2 \nabla^2 potential + # The "asymm_term" is asymmetric on x and y, so averages out to + # zero over all space + return kappa_circ * (1 - asymm_term) + (alpha_circ / grid_radii) * asymm_term - \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * - (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + @aa.grid_dec.to_array + def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + return np.zeros(shape=grid.shape[0]) - (c.f. Eliasdottir '07 eqn. A3) - In this parameterization, ra and rs are the scale radii above in angular - units (arcsec). The parameter kappa_scale is \\Sigma_0 / \\Sigma_crit. +class dPIESph(dPIE): - WARNING: This uses the "pseud-elliptical" approximation, where the ellipticity - is applied to the *potential* rather than the *mass* to ease calculation. - Use at your own risk! (And TODO Jack: fix this!) - This approximation is used by the lenstronomy PJAFFE profile (which is the - same functional form), but not by the lenstool PIEMD (also synonymous with this), - which correctly solved the differential equations for the mass-based ellipticity. - ''' def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), - ell_comps: Tuple[float, float] = (0.0, 0.0), ra: float = 0.1, rs: float = 2.0, kappa_scale: float = 0.1, ): - super(MassProfile, self).__init__(centre, ell_comps) - if ra > rs: - ra, rs = rs, ra - self.ra = ra - self.rs = rs - self.kappa_scale = kappa_scale + """ + The dual Pseudo-Isothermal Elliptical mass distribution introduced in + Eliasdottir 2007: https://arxiv.org/abs/0710.5636 - def _align_to_major_axis(self, ys, xs): - ''' - Aligns coordinates to the major axis of this halo. Returns (y', x'), - where x' is along the major axis and y' is along the minor axis. + This version is without ellipticity, so perhaps the "E" is a misnomer. - Does NOT translate, only rotates. - ''' - costheta, sintheta = self._cos_and_sin_to_x_axis() - _xs = (costheta * xs + sintheta * ys) - _ys = (-sintheta * xs + costheta * ys) - return _ys, _xs + Corresponds to a projected density profile that looks like: - def _align_from_major_axis(self, _ys, _xs): - ''' - Given _ys and _xs as offsets along the minor and major axes, - respectively, this transforms them back to the regular coordinate - system. + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) - Does NOT translate, only rotates. - ''' - costheta, sintheta = self._cos_and_sin_to_x_axis() - xs = (costheta * _xs + -sintheta * _ys) - ys = (sintheta * _xs + costheta * _ys) - return ys, xs + (c.f. Eliasdottir '07 eqn. A3) - def _ellip(self): - ellip = np.sqrt(self.ell_comps[0]**2 + self.ell_comps[1]**2) - MAX_ELLIP = 0.99999 - return min(ellip, MAX_ELLIP) + In this parameterization, ra and rs are the scale radii above in angular + units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. + + Credit: Jackson O'Donnell for writing the PyAutoLens implementation of the dPIE mass profile. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + ra + A scale radius in arc-seconds. + rs + The second scale radius in arc-seconds. + kappa_scale + Scales the overall normalization of the profile, so related to the mass + """ + + # Ensure rs > ra (things will probably break otherwise) + if ra > rs: + ra, rs = rs, ra + super().__init__(centre, (0.0, 0.0)) + self.ra = ra + self.rs = rs + self.kappa_scale = kappa_scale @aa.grid_dec.to_vector_yx + @aa.grid_dec.transform + @aa.grid_dec.relocate_to_radial_minimum def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - ys, xs = grid.T - (ycen, xcen) = self.centre - xoff, yoff = xs - xcen, ys - ycen - _ys, _xs = self._align_to_major_axis(yoff, xoff) - ellip = self._ellip() - _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) + radii = self.radial_grid_from(grid=grid, **kwargs) - # Compute the deflection magnitude of a *non-elliptical* profile - alpha_circ = self._deflection_angle(_radii) + alpha = self._deflection_angle(radii) - # This is in axes aligned to the major/minor axis - _defl_xs = alpha_circ * np.sqrt(1 - ellip) * (_xs / _radii) - _defl_ys = alpha_circ * np.sqrt(1 + ellip) * (_ys / _radii) + # now we decompose the deflection into y/x components + defl_y = alpha * grid[:,0] / radii + defl_x = alpha * grid[:,1] / radii - # And here we convert back to the real axes - defl_ys, defl_xs = self._align_from_major_axis(_defl_ys, _defl_xs) return aa.Grid2DIrregular.from_yx_1d( - defl_ys, defl_xs + defl_y, defl_x ) @aa.grid_dec.to_array + @aa.grid_dec.transform + @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - ys, xs = grid.T - (ycen, xcen) = self.centre - xoff, yoff = xs - xcen, ys - ycen - _ys, _xs = self._align_to_major_axis(yoff, xoff) - - ellip = self._ellip() - _radii = np.sqrt(_xs**2 * (1 - ellip) + _ys**2 * (1 + ellip)) - - # Compute the convergence and deflection of a *circular* profile - kappa_circ = self._convergence(_radii) - alpha_circ = self._deflection_angle(_radii) - - asymm_term = (ellip * (1 - ellip) * _xs**2 - ellip * (1 + ellip) * _ys**2) / _radii**2 - - # convergence = 1/2 \nabla \alpha = 1/2 \nabla^2 potential - # The "asymm_term" is asymmetric on x and y, so averages out to - # zero over all space - return kappa_circ * (1 - asymm_term) + (alpha_circ / _radii) * asymm_term + + # already transformed to center on profile centre so this works + radsq = (grid[:, 0]**2 + grid[:, 1]**2) + + return self._convergence(np.sqrt(radsq)) @aa.grid_dec.to_array def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): diff --git a/docs/installation/conda.rst b/docs/installation/conda.rst index 7ae7edff5..be8a89fe5 100644 --- a/docs/installation/conda.rst +++ b/docs/installation/conda.rst @@ -46,7 +46,7 @@ You may get warnings which state something like: .. code-block:: bash - ERROR: autoarray 2025.5.7.16 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible. + ERROR: autoarray 2025.5.10.1 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible. ERROR: numba 0.53.1 has requirement llvmlite<0.37,>=0.36.0rc1, but you'll have llvmlite 0.38.0 which is incompatible. If you see these messages, they do not mean that the installation has failed and the instructions below will diff --git a/docs/installation/pip.rst b/docs/installation/pip.rst index 11d1ef403..54fa54c78 100644 --- a/docs/installation/pip.rst +++ b/docs/installation/pip.rst @@ -27,7 +27,7 @@ You may get warnings which state something like: .. code-block:: bash - ERROR: autoarray 2025.5.7.16 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible. + ERROR: autoarray 2025.5.10.1 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible. ERROR: numba 0.53.1 has requirement llvmlite<0.37,>=0.36.0rc1, but you'll have llvmlite 0.38.0 which is incompatible. If you see these messages, they do not mean that the installation has failed and the instructions below will From 19e0ac58300df5d216c4c67a1b73563832e4b90b Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 13 May 2025 20:35:19 +0100 Subject: [PATCH 5/5] refactored and implemented --- .../mass/total/dual_pseudo_isothermal.py | 123 +++++++++++++----- 1 file changed, 89 insertions(+), 34 deletions(-) diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py index 949d1f9b0..3b1967f1d 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py @@ -7,27 +7,7 @@ class dPIE(MassProfile): - ''' - The dual Pseudo-Isothermal Elliptical mass distribution introduced in - Eliasdottir 2007: https://arxiv.org/abs/0710.5636 - Corresponds to a projected density profile that looks like: - - \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * - (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) - - (c.f. Eliasdottir '07 eqn. A3) - - In this parameterization, ra and rs are the scale radii above in angular - units (arcsec). The parameter kappa_scale is \\Sigma_0 / \\Sigma_crit. - - WARNING: This uses the "pseud-elliptical" approximation, where the ellipticity - is applied to the *potential* rather than the *mass* to ease calculation. - Use at your own risk! (And TODO Jack: fix this!) - This approximation is used by the lenstronomy PJAFFE profile (which is the - same functional form), but not by the lenstool PIEMD (also synonymous with this), - which correctly solved the differential equations for the mass-based ellipticity. - ''' def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -36,7 +16,42 @@ def __init__( rs: float = 2.0, kappa_scale: float = 0.1, ): - super().__init__(centre, ell_comps) + """ + The dual Pseudo-Isothermal mass profile (dPIE) without ellipticity, based on the + formulation from Eliasdottir (2007): https://arxiv.org/abs/0710.5636. + + This profile describes a circularly symmetric (non-elliptical) projected mass + distribution with two scale radii (`ra` and `rs`) and a normalization factor + `kappa_scale`. Although originally called the dPIE (Elliptical), this version + lacks ellipticity, so the "E" may be a misnomer. + + The projected surface mass density is given by: + + .. math:: + + \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * + (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) + + (See Eliasdottir 2007, Eq. A3.) + + In this implementation: + - `ra` and `rs` are scale radii in arcseconds. + - `kappa_scale` = Σ₀ / Σ_crit is the dimensionless normalization. + + Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + ra + The inner core scale radius in arcseconds. + rs + The outer truncation scale radius in arcseconds. + kappa_scale + The dimensionless normalization factor controlling the overall mass. + """ + super().__init__(centre=centre, ell_comps=ell_comps) if ra > rs: ra, rs = rs, ra @@ -81,7 +96,14 @@ def _convergence(self, radii): @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + """ + Calculate the deflection angles on a grid of (y,x) arc-second coordinates. + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the deflection angles are computed on. + """ ellip = self._ellip() grid_radii = np.sqrt(grid[:,1]**2 * (1 - ellip) + grid[:,0]**2 * (1 + ellip)) @@ -102,7 +124,17 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + """ + Returns the two dimensional projected convergence on a grid of (y,x) arc-second coordinates. + The `grid_2d_to_structure` decorator reshapes the ndarrays the convergence is outputted on. See + *aa.grid_2d_to_structure* for a description of the output. + + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the convergence is computed on. + """ ellip = self._ellip() grid_radii = np.sqrt(grid[:,1]**2 * (1 - ellip) + grid[:,0]**2 * (1 + ellip)) @@ -132,39 +164,45 @@ def __init__( kappa_scale: float = 0.1, ): """ - The dual Pseudo-Isothermal Elliptical mass distribution introduced in - Eliasdottir 2007: https://arxiv.org/abs/0710.5636 + The dual Pseudo-Isothermal mass profile (dPIE) without ellipticity, based on the + formulation from Eliasdottir (2007): https://arxiv.org/abs/0710.5636. + + This profile describes a circularly symmetric (non-elliptical) projected mass + distribution with two scale radii (`ra` and `rs`) and a normalization factor + `kappa_scale`. Although originally called the dPIE (Elliptical), this version + lacks ellipticity, so the "E" may be a misnomer. - This version is without ellipticity, so perhaps the "E" is a misnomer. + The projected surface mass density is given by: - Corresponds to a projected density profile that looks like: + .. math:: \\Sigma(R) = \\Sigma_0 (ra * rs) / (rs - ra) * (1 / \\sqrt(ra^2 + R^2) - 1 / \\sqrt(rs^2 + R^2)) - (c.f. Eliasdottir '07 eqn. A3) + (See Eliasdottir 2007, Eq. A3.) - In this parameterization, ra and rs are the scale radii above in angular - units (arcsec). The parameter `kappa_scale` is \\Sigma_0 / \\Sigma_crit. + In this implementation: + - `ra` and `rs` are scale radii in arcseconds. + - `kappa_scale` = Σ₀ / Σ_crit is the dimensionless normalization. - Credit: Jackson O'Donnell for writing the PyAutoLens implementation of the dPIE mass profile. + Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. ra - A scale radius in arc-seconds. + The inner core scale radius in arcseconds. rs - The second scale radius in arc-seconds. + The outer truncation scale radius in arcseconds. kappa_scale - Scales the overall normalization of the profile, so related to the mass + The dimensionless normalization factor controlling the overall mass. """ # Ensure rs > ra (things will probably break otherwise) if ra > rs: ra, rs = rs, ra - super().__init__(centre, (0.0, 0.0)) + super().__init__(centre=centre, ell_comps=(0.0, 0.0)) self.ra = ra self.rs = rs self.kappa_scale = kappa_scale @@ -173,7 +211,14 @@ def __init__( @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): + """ + Calculate the deflection angles on a grid of (y,x) arc-second coordinates. + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the deflection angles are computed on. + """ radii = self.radial_grid_from(grid=grid, **kwargs) alpha = self._deflection_angle(radii) @@ -190,7 +235,17 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - + """ + Returns the two dimensional projected convergence on a grid of (y,x) arc-second coordinates. + + The `grid_2d_to_structure` decorator reshapes the ndarrays the convergence is outputted on. See + *aa.grid_2d_to_structure* for a description of the output. + + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the convergence is computed on. + """ # already transformed to center on profile centre so this works radsq = (grid[:, 0]**2 + grid[:, 1]**2)