From 025ac7ac02092081618b6f889e9349866d748870 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 09:56:35 +0100 Subject: [PATCH 1/9] remove relocate_to_radial_minimum --- autogalaxy/profiles/mass/total/pseudo_isothermal.py | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/autogalaxy/profiles/mass/total/pseudo_isothermal.py b/autogalaxy/profiles/mass/total/pseudo_isothermal.py index 877ad5eeb..d90587cd7 100644 --- a/autogalaxy/profiles/mass/total/pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/pseudo_isothermal.py @@ -254,7 +254,6 @@ def _ellip(self): @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): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. @@ -278,7 +277,6 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): ) @aa.grid_dec.transform - @aa.grid_dec.relocate_to_radial_minimum def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. @@ -369,7 +367,6 @@ def _ellip(self): @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): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. @@ -395,7 +392,6 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): ) @aa.grid_dec.transform - @aa.grid_dec.relocate_to_radial_minimum def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. @@ -489,7 +485,6 @@ def __init__( @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): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. @@ -526,8 +521,8 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): return self.rotated_grid_from_reference_frame_from( grid=np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T), **kwargs ) + @aa.grid_dec.transform - @aa.grid_dec.relocate_to_radial_minimum def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. @@ -689,7 +684,6 @@ 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): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. @@ -718,7 +712,6 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @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): """ Returns the two dimensional projected convergence on a grid of (y,x) arc-second coordinates. @@ -754,7 +747,6 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): return np.zeros(shape=grid.shape[0]) - class dPIEPSph(dPIEP): def __init__( self, @@ -814,7 +806,6 @@ def __init__( @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): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. @@ -836,7 +827,6 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @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): """ Returns the two dimensional projected convergence on a grid of (y,x) arc-second coordinates. From 093b72ca680fec6330d9773c5e12f1814f52746a Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:11:24 +0100 Subject: [PATCH 2/9] update dPIEP unit tests and move to separate module --- autogalaxy/operate/image.py | 1 + autogalaxy/profiles/mass/total/__init__.py | 4 +- ...rmal.py => dual_pseudo_isothermal_mass.py} | 444 ++++-------------- ...py => dual_pseudo_isothermal_potential.py} | 68 +-- ... test_dual_pseudo_isothermal_potential.py} | 38 +- 5 files changed, 158 insertions(+), 397 deletions(-) rename autogalaxy/profiles/mass/total/{pseudo_isothermal.py => dual_pseudo_isothermal_mass.py} (62%) rename autogalaxy/profiles/mass/total/{dual_pseudo_isothermal.py => dual_pseudo_isothermal_potential.py} (75%) rename test_autogalaxy/profiles/mass/total/{test_dual_pseudo_isothermal.py => test_dual_pseudo_isothermal_potential.py} (62%) diff --git a/autogalaxy/operate/image.py b/autogalaxy/operate/image.py index 53aadcdff..89994814d 100644 --- a/autogalaxy/operate/image.py +++ b/autogalaxy/operate/image.py @@ -195,6 +195,7 @@ def visibilities_from( return aa.Visibilities.zeros(shape_slim=(transformer.uv_wavelengths.shape[0],)) + class OperateImageList(OperateImage): """ Packages methods which operate on the list of 2D images returned from the `image_2d_list_from` function of a light diff --git a/autogalaxy/profiles/mass/total/__init__.py b/autogalaxy/profiles/mass/total/__init__.py index 6e4b00185..bcbde48da 100644 --- a/autogalaxy/profiles/mass/total/__init__.py +++ b/autogalaxy/profiles/mass/total/__init__.py @@ -1,5 +1,5 @@ -from .dual_pseudo_isothermal import dPIE, dPIESph -from .pseudo_isothermal import PIEMD, dPIEMD, dPIEMDSph, dPIEP, dPIEPSph +from .dual_pseudo_isothermal_potential import dPIEP, dPIEPSph +from .dual_pseudo_isothermal_mass import PIEMD, dPIEMD, dPIEMDSph from .isothermal import Isothermal, IsothermalSph from .isothermal_core import IsothermalCore, IsothermalCoreSph from .power_law import PowerLaw, PowerLawSph diff --git a/autogalaxy/profiles/mass/total/pseudo_isothermal.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py similarity index 62% rename from autogalaxy/profiles/mass/total/pseudo_isothermal.py rename to autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index d90587cd7..853ce1442 100644 --- a/autogalaxy/profiles/mass/total/pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -5,11 +5,12 @@ from autogalaxy.profiles.mass.abstract.abstract import MassProfile # Within this profile family, PIEMD, dPIEMD, and dPIEMDSph are directly ported from Lenstool's C code, and have been thoroughly annotated and adapted for PyAutoLens. -# The dPIEP and dPIEPSph profiles are modified from the original `dPIE` and `dPIESph`, which were implemented to PyAutoLens by Jackson O'Donnell. +# The dPIEP and dPIEPSph profiles are modified from the original `dPIEP` and `dPIEPSph`, which were implemented to PyAutoLens by Jackson O'Donnell. + def _ci05(x, y, eps, rcore): """ - Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / E0 for PIEMD at any positions (x,y), + Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / E0 for PIEMD at any positions (x,y), see Kassiola & Kovner(1993) Eq. 4.1.2, which is the integral of Eq. 2.3.8. Note here b0(or called E0) is out of the `_ci05`. @@ -35,31 +36,35 @@ def _ci05(x, y, eps, rcore): # Define intermediate complex variables zci = np.complex128(complex(0.0, -0.5 * (1.0 - eps * eps) / sqe)) - znum = np.complex128(axis_ratio * x + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio)) + znum = np.complex128( + axis_ratio * x + + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio) + ) zden = np.complex128(x + 1j * (2.0 * rcore * sqe - y)) - + # zis = znum / zden = (a+bi)/(c+di) = [(ac+bd)+(bc-ad i)] / (c^2+d^2) norm = zden.real * zden.real + zden.imag * zden.imag # |zden|^2 zis_re = (znum.real * zden.real + znum.imag * zden.imag) / norm zis_im = (znum.imag * zden.real - znum.real * zden.imag) / norm zis = np.complex128(zis_re + 1j * zis_im) - + # ln(zis) = ln(|zis|) + i*Arg(zis) zis_mag = np.abs(zis) zis_re = np.log(zis_mag) zis_im = np.angle(zis) zis = np.complex128(zis_re + 1j * zis_im) - + # I'* = zres = zci * ln(zis) zres = zci * zis - + return zres + def _ci05f(x, y, eps, rcore, rcut): """ - Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / (b0 * ra / (rs - ra)) for dPIEMD at any positions (x,y), - which is the integral of Eq. 2.3.8 in Kassiola & Kovner(1993). - + Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / (b0 * ra / (rs - ra)) for dPIEMD at any positions (x,y), + which is the integral of Eq. 2.3.8 in Kassiola & Kovner(1993). + Note here (b0 * ra / (rs - ra)) is out of the `_ci05f`. The only difference of integral of Eq. 2.3.8 between dPIEMD and PIEMD is the \\kappa: \\kappa(r_{em})_{dPIEMD} = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}). I*_{dPIEMD} = ra / (rs - ra) * (I*_{PIEMD}(ra) - I*_{PIEMD}(ra)) @@ -89,10 +94,15 @@ def _ci05f(x, y, eps, rcore, rcut): # Define intermediate complex variables zci = np.complex128(complex(0.0, -0.5 * (1.0 - eps * eps) / sqe)) - znum_rc = np.complex128(axis_ratio * x + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio)) # a + bi - zden_rc = np.complex128(x + 1j * (2.0 * rcore * sqe - y)) # c + di - znum_rcut = np.complex128(axis_ratio * x + 1j * (2.0 * sqe * np.sqrt(rcut * rcut + rem2) - y / axis_ratio)) # a + ei - zden_rcut = np.complex128(x + 1j * (2.0 * rcut * sqe - y)) # c + fi + znum_rc = np.complex128( + axis_ratio * x + + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio) + ) # a + bi + zden_rc = np.complex128(x + 1j * (2.0 * rcore * sqe - y)) # c + di + znum_rcut = np.complex128( + axis_ratio * x + 1j * (2.0 * sqe * np.sqrt(rcut * rcut + rem2) - y / axis_ratio) + ) # a + ei + zden_rcut = np.complex128(x + 1j * (2.0 * rcut * sqe - y)) # c + fi # zis_rc = znum_rc / zden_rc = (a+bi)/(c+di) # zis_rcut = znum_rcut / zden_rcut = (a+ei)/(c+fi) @@ -102,10 +112,10 @@ def _ci05f(x, y, eps, rcore, rcut): # = (aa + bb*i) * (cc -dd*i) / (cc^2 + dd^2) # = [(aa*cc + bb*dd) / (cc^ + dd^2)] + [(bb*cc - aa*dd) / (cc^2 + dd^2)]*i # = aaa + bbb*i - aa = znum_rc.real * zden_rc.real - znum_rc.imag * zden_rcut.imag # ac - bf - bb = znum_rc.real * zden_rcut.imag + znum_rc.imag * zden_rc.real # af + bc - cc = znum_rc.real * zden_rc.real - zden_rc.imag * znum_rcut.imag # ac - de - dd = znum_rc.real * zden_rc.imag + zden_rc.real * znum_rcut.imag # ad + ce + aa = znum_rc.real * zden_rc.real - znum_rc.imag * zden_rcut.imag # ac - bf + bb = znum_rc.real * zden_rcut.imag + znum_rc.imag * zden_rc.real # af + bc + cc = znum_rc.real * zden_rc.real - zden_rc.imag * znum_rcut.imag # ac - de + dd = znum_rc.real * zden_rc.imag + zden_rc.real * znum_rcut.imag # ad + ce norm = cc * cc + dd * dd aaa = (aa * cc + bb * dd) / norm bbb = (bb * cc - aa * dd) / norm @@ -122,6 +132,7 @@ def _ci05f(x, y, eps, rcore, rcut): return zres + def _mdci05(x, y, eps, rcore, b0): """ Returns the second derivatives (Hessian matrix) of the lens potential as complex number for PIEMD at any positions (x,y): @@ -152,13 +163,13 @@ def _mdci05(x, y, eps, rcore, b0): cxro = (1.0 + eps) * (1.0 + eps) cyro = (1.0 - eps) * (1.0 - eps) ci = 0.5 * (1.0 - eps * eps) / sqe - wrem = np.sqrt(rcore * rcore + x * x / cxro + y * y / cyro) # √(w(x,y)) + wrem = np.sqrt(rcore * rcore + x * x / cxro + y * y / cyro) # √(w(x,y)) num1 = 2.0 * sqe * wrem - y * axis_ratio_inv - den1 = axis_ratio * axis_ratio * x * x + num1 * num1 # |q * x + num1*i|^2 + den1 = axis_ratio * axis_ratio * x * x + num1 * num1 # |q * x + num1*i|^2 num2 = 2.0 * rcore * sqe - y - den2 = x * x + num2 * num2 # |x + num2*i|^2 + den2 = x * x + num2 * num2 # |x + num2*i|^2 - # eg. + # eg. # ∂²ψ/∂x² = Re(∂I*/∂x) = b0 * didxre # ∂I*/∂x = b0 * ci * (-i) * ∂(ln{u(x,y)} - ln{v(x,y)})∂x # = b0 * ci * (-i) * (1/u * ∂u/∂x - 1/v * ∂v/∂x) @@ -171,7 +182,7 @@ def _mdci05(x, y, eps, rcore, b0): # = {q + [2.0 * sqe * x / cxro / wrem] * i} * {q * x - (2.0 * sqe * wrem - y / q)*i} / den1 # = {q^2 * x + 4.0 * sqe^2 * x - y / q * 2.0 * sqe * x / cxro / wrem} / den1 + q * { (2.0 * sqe * x^2 / cxro / wrem) - (2.0 * sqe * wrem - y / q)} / den1 * i # = {x - 2.0 * sqe * x * y * q / cyro / wrem} / den1 + q * { (2.0 * sqe * x^2 / cxro / wrem) - (2.0 * sqe * wrem - y / q)} / den1 * i - # (-i) * (1/u * ∂u/∂x) = (2.0 * sqe * x * y * q / cyro / wrem - x) / den1 * i + # (-i) * (1/u * ∂u/∂x) = (2.0 * sqe * x * y * q / cyro / wrem - x) / den1 * i # + q * { (2.0 * sqe * x^2 / cxro / wrem) - (2.0 * sqe * wrem - y / q)} / den1 # ∂v/∂x = 1 + ∂(num2)/∂x * i # = 1 @@ -184,16 +195,20 @@ def _mdci05(x, y, eps, rcore, b0): # Compute second derivatives didxre = ci * ( - axis_ratio * (2.0 * sqe * x * x / cxro / wrem - 2.0 * sqe * wrem + y * axis_ratio_inv) / den1 + axis_ratio + * (2.0 * sqe * x * x / cxro / wrem - 2.0 * sqe * wrem + y * axis_ratio_inv) + / den1 + num2 / den2 ) - didyre = ci * ( - (2.0 * sqe * x * y * axis_ratio / cyro / wrem - x) / den1 - + x / den2 - ) + didyre = ci * ((2.0 * sqe * x * y * axis_ratio / cyro / wrem - x) / den1 + x / den2) didyim = ci * ( - (2.0 * sqe * wrem * axis_ratio_inv - y * axis_ratio_inv * axis_ratio_inv - 4.0 * eps * y / cyro - + 2.0 * sqe * y * y / cyro / wrem * axis_ratio_inv) / den1 + ( + 2.0 * sqe * wrem * axis_ratio_inv + - y * axis_ratio_inv * axis_ratio_inv + - 4.0 * eps * y / cyro + + 2.0 * sqe * y * y / cyro / wrem * axis_ratio_inv + ) + / den1 - num2 / den2 ) @@ -203,7 +218,8 @@ def _mdci05(x, y, eps, rcore, b0): c = b0 * didyre # ∂²ψ/∂y∂x d = b0 * didyim # ∂²ψ/∂y² - return a,b,c,d + return a, b, c, d + class PIEMD(MassProfile): def __init__( @@ -214,22 +230,22 @@ def __init__( b0: float = 0.1, ): """ - The Pseudo Isothermal Elliptical Mass Distribution(PIEMD) profiles, based on the formulaiton from + The Pseudo Isothermal Elliptical Mass Distribution(PIEMD) profiles, based on the formulaiton from Kassiola & Kovner(1993) https://articles.adsabs.harvard.edu/pdf/1993ApJ...417..450K. This profile is ported from Lenstool's C code, which has the same formulation. - This proflie describes an elliptic isothermal mass distribution with a finite core: + This proflie describes an elliptic isothermal mass distribution with a finite core: \\rho \\propto (ra^2 + R^2)^{-1} The convergence is given by: \\kappa(r_{em}) = \\kappa_0 * ra / \\sqrt{ ra^2 + r_{em}^2 } (see Kassiola & Kovner(1993), Eq. 4.1.1) where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2, (see Kassiola & Kovner(1993), Eq. 2.3.6) - and \\kappa_0 = b_0 / 2 / r_a. + and \\kappa_0 = b_0 / 2 / r_a. In this implementation: - `ra` is the core radius in unit of arcseconds. - - `b0` is the lens strength in unit of arcseconds, when ra->0 & q->1, b0 is the Einstein radius. + - `b0` is the lens strength in unit of arcseconds, when ra->0 & q->1, b0 is the Einstein radius. `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}). `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 @@ -246,12 +262,12 @@ def __init__( self.ra = ra self.b0 = b0 - + 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 @aa.grid_dec.transform def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @@ -273,11 +289,12 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # And here we convert back to the real axes return self.rotated_grid_from_reference_frame_from( - grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), **kwargs + grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), + **kwargs, ) @aa.grid_dec.transform - def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): + def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. @@ -296,17 +313,18 @@ def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): ) return hessian_yy, hessian_xy, hessian_yx, hessian_xx - - def analytical_magnification_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): - hessian_yy, hessian_xy, hessian_yx, hessian_xx = self.analytical_hessian_2d_from( - grid=grid + def analytical_magnification_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): + + hessian_yy, hessian_xy, hessian_yx, hessian_xx = ( + self.analytical_hessian_2d_from(grid=grid) ) det_A = (1 - hessian_xx) * (1 - hessian_yy) - hessian_xy * hessian_yx return aa.Array2D(values=1.0 / det_A, mask=grid.mask) - + + class dPIEMD(MassProfile): def __init__( self, @@ -317,26 +335,26 @@ def __init__( b0: float = 0.1, ): """ - The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles, which is a *two component PIEMD* with both a core radius and a truncation radius, + The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles, which is a *two component PIEMD* with both a core radius and a truncation radius, see Eliasdottir (2007): https://arxiv.org/abs/0710.5636 This profile is ported from Lenstool's C code, which has the same formulation. - This proflie describes an elliptic isothermal mass distribution with a finite core, \\rho \~ r^{-2} while in the transition region (ra<=R<=rs), - and \\rho \~ r^{-4} in the outer parts: + This proflie describes an elliptic isothermal mass distribution with a finite core, \\rho \~ r^{-2} while in the transition region (ra<=R<=rs), + and \\rho \~ r^{-4} in the outer parts: \\rho \\propto [(ra^2 + R^2) (rs^2 + R^2)]^{-1} - The convergence is given by two PIEMD with core radius ra and rs: + The convergence is given by two PIEMD with core radius ra and rs: \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) - where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. - Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIE}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. - There is \\frac{\\sigma_{dPIE}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, + where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. + Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEP}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. + There is \\frac{\\sigma_{dPIEP}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, thus E0(Kassiola & Kovner(1993)) = b0 = E0(Eliasdottir (2007)) * (rs^2 - ra^2) / rs^2. So when s->\\infty and a->0, they are equivalent. In this implementation: - `ra` is the core radius in unit of arcseconds. - `rs` is the truncation radius in unit of arcseconds. - - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. + - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 @@ -359,12 +377,12 @@ def __init__( self.ra = ra self.rs = rs self.b0 = b0 - + 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 @aa.grid_dec.transform def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): @@ -378,9 +396,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): """ ellip = self._ellip() factor = self.b0 * self.rs / (self.rs - self.ra) - zis = _ci05f( - x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, rcut=self.rs - ) + zis = _ci05f(x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, rcut=self.rs) # This is in axes aligned to the major/minor axis deflection_x = zis.real @@ -388,11 +404,12 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # And here we convert back to the real axes return self.rotated_grid_from_reference_frame_from( - grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), **kwargs + grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), + **kwargs, ) - + @aa.grid_dec.transform - def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): + def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. Hessian_dPIEMD = rs * (rs - ra) * ( Hessian_PIEMD(ra) - Hessian_PIEMD(rs)) @@ -406,7 +423,7 @@ def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): if grid.ndim != 2 or grid.shape[1] != 2: raise ValueError("Grid must be a 2D array with shape (n, 2)") ellip = self._ellip() - + t05 = self.rs / (self.rs - self.ra) g05c_a, g05c_b, g05c_c, g05c_d = _mdci05( x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, b0=self.b0 @@ -422,46 +439,47 @@ def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): hessian_yy = t05 * (g05c_d - g05cut_d) return hessian_yy, hessian_xy, hessian_yx, hessian_xx - - def analytical_magnification_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): - - hessian_yy, hessian_xy, hessian_yx, hessian_xx = self.analytical_hessian_2d_from( - grid=grid + + def analytical_magnification_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): + + hessian_yy, hessian_xy, hessian_yx, hessian_xx = ( + self.analytical_hessian_2d_from(grid=grid) ) det_A = (1 - hessian_xx) * (1 - hessian_yy) - hessian_xy * hessian_yx return aa.Array2D(values=1.0 / det_A, mask=grid.mask) - + + class dPIEMDSph(dPIEMD): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), ra: float = 0.1, rs: float = 2.0, - b0: float = 1.0 + b0: float = 1.0, ): """ - The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles without ellipticity, which is a *two component PIEMD* with both a core radius and a truncation radius, + The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles without ellipticity, which is a *two component PIEMD* with both a core radius and a truncation radius, see Eliasdottir (2007): https://arxiv.org/abs/0710.5636 This profile is ported from Lenstool's C code, which has the same formulation. - This proflie describes an spherical isothermal mass distribution with a finite core, \\rho \~ r^{-2} while in the transition region (ra<=R<=rs), - and \\rho \~ r^{-4} in the outer parts: + This proflie describes an spherical isothermal mass distribution with a finite core, \\rho \~ r^{-2} while in the transition region (ra<=R<=rs), + and \\rho \~ r^{-4} in the outer parts: \\rho \\propto [(ra^2 + R^2) (rs^2 + R^2)]^{-1} - The convergence is given by two PIEMD with core radius ra and rs: + The convergence is given by two PIEMD with core radius ra and rs: \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) - where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. - Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIE}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. - There is \\frac{\\sigma_{dPIE}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, + where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. + Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEP}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. + There is \\frac{\\sigma_{dPIEP}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, thus E0(Kassiola & Kovner(1993)) = b0 = E0(Eliasdottir (2007)) * (rs^2 - ra^2) / rs^2. So when s->\\infty and a->0, they are equivalent. In this implementation: - `ra` is the core radius in unit of arcseconds. - `rs` is the truncation radius in unit of arcseconds. - - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. + - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 @@ -488,7 +506,7 @@ def __init__( def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): """ Calculate the deflection angles on a grid of (y,x) arc-second coordinates. - Faster and equivalent to Eliasdottir (2007), see Eq. A19 and Eq. A20. + Faster and equivalent to Eliasdottir (2007), see Eq. A19 and Eq. A20. f(R,a,s) = {R/a} / {1 + \\sqrt{1 + (R/a)^2}} - {R/s} / {1 + \\sqrt{1 + (R/s)^2}} = R / {\\sqrt{a^2 + R^2} + a} - R / {\\sqrt{s^2 + R^2} + s} @@ -523,10 +541,10 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): ) @aa.grid_dec.transform - def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): + def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. - Chain rule of second derivatives: + Chain rule of second derivatives: ∂²ψ/∂x² = ∂²ψ/∂R² * (∂R/∂x)² + ∂²R/∂x² * ∂ψ/∂R ∂²ψ/∂y² = ∂²ψ/∂R² * (∂R/∂y)² + ∂²R/∂y² * ∂ψ/∂R ∂²ψ/∂x∂y = ∂²ψ/∂R² * ∂R/∂x * ∂R/∂y + ∂²R/∂x∂y * ∂ψ/∂R @@ -564,10 +582,12 @@ def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): # = p R2 = grid[:, 1] * grid[:, 1] + grid[:, 0] * grid[:, 0] z = np.sqrt(R2 + a * a) - a - np.sqrt(R2 + s * s) + s - p = (1.0 - a / np.sqrt(a * a + R2)) * a / R2 - (1.0 - s / np.sqrt(s * s + R2)) * s / R2 - X = grid[:, 1] * grid[:, 1] / R2 # x^2 / R^2 - Y = grid[:, 0] * grid[:, 0] / R2 # y^2 / R^2 - XY = grid[:, 1] * grid[:, 0] / R2 # x*y / R^2 + p = (1.0 - a / np.sqrt(a * a + R2)) * a / R2 - ( + 1.0 - s / np.sqrt(s * s + R2) + ) * s / R2 + X = grid[:, 1] * grid[:, 1] / R2 # x^2 / R^2 + Y = grid[:, 0] * grid[:, 0] / R2 # y^2 / R^2 + XY = grid[:, 1] * grid[:, 0] / R2 # x*y / R^2 # ∂²ψ/∂x² = ∂²ψ/∂R² * (∂R/∂x)² + ∂²R/∂x² * ∂ψ/∂R # = p * (x / R)^2 + y^2 / R^3 * z / R @@ -589,261 +609,3 @@ def analytical_hessian_2d_from(self, grid: 'aa.type.Grid2DLike', **kwargs): hessian_yy = t05 * (p * Y + z * X / R2) return hessian_yy, hessian_xy, hessian_yx, hessian_xx - -class dPIEP(MassProfile): - - 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, - b0: float = 1.0, - ): - """ - The dual Pseudo Isothermal Elliptical Potential (dPIEP) with pseudo-ellipticity on potential, 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` is the core radius in unit of arcseconds. - - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. - `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) - `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 - - Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. - Note: To ensure consistency, kappa_scale was replaced with b0, and the corresponding code was adjusted accordingly. - - 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. - b0 - The lens strength in arcseconds. - """ - super().__init__(centre=centre, ell_comps=ell_comps) - - if ra > rs: - ra, rs = rs, ra - - self.ra = ra - self.rs = rs - self.b0 = b0 - - 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. - """ - a, s = self.ra, self.rs - radii = np.maximum(radii, 1e-8) - f = ( - radii / (a + np.sqrt(a**2 + radii**2)) - - radii / (s + np.sqrt(s**2 + radii**2)) - ) - - # c.f. Eliasdottir '07 eq. A23 - # magnitude of deflection - # alpha = self.E0 * (s + a) / s * f - alpha = self.b0 * s / (s - a) * f - return alpha - - def _convergence(self, radii): - radsq = radii * radii - a, s = self.ra, self.rs - # c.f. Eliasdottir '07 eqn (A3) - # return ( - # self.E0 / 2 * (s + a) / s * - # (1/np.sqrt(a**2 + radsq) - 1/np.sqrt(s**2 + radsq)) - # ) - return ( - self.b0 / 2 * s / (s - a) * - (1/np.sqrt(a**2 + radsq) - 1/np.sqrt(s**2 + radsq)) - ) - - @aa.grid_dec.to_vector_yx - @aa.grid_dec.transform - 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) - ) - - # 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_vector_yx - @aa.grid_dec.transform - 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) - ) - - # Compute the convergence and deflection of a *circular* profile - kappa_circ = self._convergence(grid_radii) - alpha_circ = self._deflection_angle(grid_radii) - - asymm_term = ( - ellip * (1 - ellip) * grid[:, 1] ** 2 - - ellip * (1 + ellip) * grid[:, 0] ** 2 - ) / grid_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 / grid_radii) * asymm_term - - @aa.grid_dec.to_array - def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): - return np.zeros(shape=grid.shape[0]) - -class dPIEPSph(dPIEP): - def __init__( - self, - centre: Tuple[float, float] = (0.0, 0.0), - ra: float = 0.1, - rs: float = 2.0, - b0: float = 1.0, - ): - """ - 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` is the core radius in unit of arcseconds. - - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. - `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) - `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 - - Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. - Note: This dPIEPSph should be the same with dPIEMDSph for their same mathamatical formulations. - - 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. - b0 - The lens strength in arcseconds. - """ - - # Ensure rs > ra (things will probably break otherwise) - if ra > rs: - ra, rs = rs, ra - - super().__init__(centre=centre, ell_comps=(0.0, 0.0)) - - self.ra = ra - self.rs = rs - self.b0 = b0 - - @aa.grid_dec.to_vector_yx - @aa.grid_dec.transform - 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) - - # now we decompose the deflection into y/x components - defl_y = alpha * grid[:, 0] / radii - defl_x = alpha * grid[:, 1] / radii - - return aa.Grid2DIrregular.from_yx_1d(defl_y, defl_x) - - @aa.grid_dec.to_array - @aa.grid_dec.transform - 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 - - 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]) diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py similarity index 75% rename from autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py rename to autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py index dd29d739d..99d637963 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py @@ -5,7 +5,7 @@ from autogalaxy.profiles.mass.abstract.abstract import MassProfile -class dPIE(MassProfile): +class dPIEP(MassProfile): def __init__( self, @@ -13,15 +13,15 @@ def __init__( ell_comps: Tuple[float, float] = (0.0, 0.0), ra: float = 0.1, rs: float = 2.0, - kappa_scale: float = 0.1, + b0: float = 1.0, ): """ - The dual Pseudo-Isothermal mass profile (dPIE) without ellipticity, based on the + The dual Pseudo Isothermal Elliptical Potential (dPIEP) with pseudo-ellipticity on potential, 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 + `kappa_scale`. Although originally called the dPIEP (Elliptical), this version lacks ellipticity, so the "E" may be a misnomer. The projected surface mass density is given by: @@ -34,10 +34,13 @@ def __init__( (See Eliasdottir 2007, Eq. A3.) In this implementation: - - `ra` and `rs` are scale radii in arcseconds. - - `kappa_scale` = Σ₀ / Σ_crit is the dimensionless normalization. + - `ra` is the core radius in unit of arcseconds. + - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. + `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) + `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. + Note: To ensure consistency, kappa_scale was replaced with b0, and the corresponding code was adjusted accordingly. Parameters ---------- @@ -47,8 +50,8 @@ def __init__( 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. + b0 + The lens strength in arcseconds. """ super().__init__(centre=centre, ell_comps=ell_comps) @@ -57,7 +60,7 @@ def __init__( self.ra = ra self.rs = rs - self.kappa_scale = kappa_scale + self.b0 = b0 def _ellip(self): ellip = jnp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) @@ -66,28 +69,29 @@ def _ellip(self): def _deflection_angle(self, radii): """ - For a circularly symmetric dPIE profile, computes the magnitude of the deflection at each radius. + For a circularly symmetric dPIEP 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 + jnp.sqrt(1 + r_ra * r_ra)) - r_rs / ( - 1 + jnp.sqrt(1 + r_rs * r_rs) + a, s = self.ra, self.rs + radii = jnp.maximum(radii, 1e-8) + f = radii / (a + jnp.sqrt(a**2 + radii**2)) - radii / ( + s + jnp.sqrt(s**2 + radii**2) ) - ra, rs = self.ra, self.rs - # c.f. Eliasdottir '07 eq. A19 + # c.f. Eliasdottir '07 eq. A23 # magnitude of deflection - alpha = 2 * self.kappa_scale * ra * rs / (rs - ra) * f + # alpha = self.E0 * (s + a) / s * f + alpha = self.b0 * s / (s - a) * 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) + self.b0 + / 2 + * s / (s - a) * (1 / jnp.sqrt(a**2 + radsq) - 1 / jnp.sqrt(s**2 + radsq)) ) @@ -162,22 +166,21 @@ def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): return jnp.zeros(shape=grid.shape[0]) -class dPIESph(dPIE): - +class dPIEPSph(dPIEP): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), ra: float = 0.1, rs: float = 2.0, - kappa_scale: float = 0.1, + b0: float = 1.0, ): """ - The dual Pseudo-Isothermal mass profile (dPIE) without ellipticity, based on the + The dual Pseudo-Isothermal mass profile (dPIEP) 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 + `kappa_scale`. Although originally called the dPIEP (Elliptical), this version lacks ellipticity, so the "E" may be a misnomer. The projected surface mass density is given by: @@ -190,10 +193,13 @@ def __init__( (See Eliasdottir 2007, Eq. A3.) In this implementation: - - `ra` and `rs` are scale radii in arcseconds. - - `kappa_scale` = Σ₀ / Σ_crit is the dimensionless normalization. + - `ra` is the core radius in unit of arcseconds. + - `b0` is the lens strength in unit of arcseconds, when ra->0 & rs->\\infty & q->1, b0 is the Einstein radius. + `b0` is related to the central velocity dispersion \\sigma_0: b_0 = 4\\pi * \\sigma_0^2 / c^2 * (D_{LS} / D_{S}) + `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. + Note: This dPIEPSph should be the same with dPIEMDSph for their same mathamatical formulations. Parameters ---------- @@ -203,8 +209,8 @@ def __init__( 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. + b0 + The lens strength in arcseconds. """ # Ensure rs > ra (things will probably break otherwise) @@ -215,7 +221,7 @@ def __init__( self.ra = ra self.rs = rs - self.kappa_scale = kappa_scale + self.b0 = b0 @aa.grid_dec.to_vector_yx @aa.grid_dec.transform diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py similarity index 62% rename from test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py rename to test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py index 0bc7d06ed..dfe23df84 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py @@ -6,42 +6,38 @@ 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) + mp = ag.mp.dPIEPSph(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) - mp = ag.mp.dPIESph(centre=(-0.1, 0.1), kappa_scale=5.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPSph(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) - mp = ag.mp.dPIE( - centre=(0, 0), ell_comps=(0.0, 0.333333), kappa_scale=1.0, ra=2.0, rs=3.0 - ) + mp = ag.mp.dPIEP(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.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 - ) + mp = ag.mp.dPIEP(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.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 + elliptical = ag.mp.dPIEP( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.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) + spherical = ag.mp.dPIEPSph(centre=(1.1, 1.1), b0=12.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 @@ -52,40 +48,36 @@ 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) + mp = ag.mp.dPIEPSph(centre=(0.0, 0.0), b0=8.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 - ) + mp = ag.mp.dPIEP(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=4.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 - ) + mp = ag.mp.dPIEP(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=8.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 + mp = ag.mp.dPIEP( + centre=(0.0, 0.0), ell_comps=(0.0, 0.333333), b0=4.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 + elliptical = ag.mp.dPIEP( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.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) + spherical = ag.mp.dPIEPSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) assert elliptical.convergence_2d_from(grid=grid).array == pytest.approx( spherical.convergence_2d_from(grid=grid).array, 1e-4 From 2b866fe0ec684d3d78ed3dd47d46a0c2baf1de35 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:13:55 +0100 Subject: [PATCH 3/9] add prior config file --- .../dual_psuedo_isothermal_potential.yaml | 126 ++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml diff --git a/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml b/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml new file mode 100644 index 000000000..4fe6d9e84 --- /dev/null +++ b/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml @@ -0,0 +1,126 @@ +dPIEP: + 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 +dPIEPSph: + 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 From 4635e7e2ded2397f56e10c7ff7a792bfc741fe56 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:14:38 +0100 Subject: [PATCH 4/9] dPIEP -> dPIEPotential --- autogalaxy/profiles/mass/__init__.py | 4 ++-- autogalaxy/profiles/mass/total/__init__.py | 2 +- .../mass/total/dual_pseudo_isothermal_mass.py | 10 ++++---- .../total/dual_pseudo_isothermal_potential.py | 16 ++++++------- .../test_dual_pseudo_isothermal_potential.py | 24 +++++++++---------- 5 files changed, 28 insertions(+), 28 deletions(-) diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index c13e2d6c6..5aaa45e70 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -4,8 +4,8 @@ dPIEMD, dPIEMDSph, PIEMD, - dPIEP, - dPIEPSph, + dPIEPotential, + dPIEPotentialSph, PowerLawCore, PowerLawCoreSph, PowerLawBroken, diff --git a/autogalaxy/profiles/mass/total/__init__.py b/autogalaxy/profiles/mass/total/__init__.py index bcbde48da..1b0ecc88b 100644 --- a/autogalaxy/profiles/mass/total/__init__.py +++ b/autogalaxy/profiles/mass/total/__init__.py @@ -1,4 +1,4 @@ -from .dual_pseudo_isothermal_potential import dPIEP, dPIEPSph +from .dual_pseudo_isothermal_potential import dPIEPotential, dPIEPotentialSph from .dual_pseudo_isothermal_mass import PIEMD, dPIEMD, dPIEMDSph from .isothermal import Isothermal, IsothermalSph from .isothermal_core import IsothermalCore, IsothermalCoreSph diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index 853ce1442..be2ccd09f 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -5,7 +5,7 @@ from autogalaxy.profiles.mass.abstract.abstract import MassProfile # Within this profile family, PIEMD, dPIEMD, and dPIEMDSph are directly ported from Lenstool's C code, and have been thoroughly annotated and adapted for PyAutoLens. -# The dPIEP and dPIEPSph profiles are modified from the original `dPIEP` and `dPIEPSph`, which were implemented to PyAutoLens by Jackson O'Donnell. +# The dPIEPotential and dPIEPotentialSph profiles are modified from the original `dPIEPotential` and `dPIEPotentialSph`, which were implemented to PyAutoLens by Jackson O'Donnell. def _ci05(x, y, eps, rcore): @@ -347,8 +347,8 @@ def __init__( \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. - Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEP}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. - There is \\frac{\\sigma_{dPIEP}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, + Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEPotential}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. + There is \\frac{\\sigma_{dPIEPotential}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, thus E0(Kassiola & Kovner(1993)) = b0 = E0(Eliasdottir (2007)) * (rs^2 - ra^2) / rs^2. So when s->\\infty and a->0, they are equivalent. In this implementation: @@ -472,8 +472,8 @@ def __init__( \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. - Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEP}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. - There is \\frac{\\sigma_{dPIEP}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, + Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEPotential}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. + There is \\frac{\\sigma_{dPIEPotential}^2}{\\sigma_0^2} = \\frac{2}{3} \frac{rs^2}{rs^2-ra^2}, thus E0(Kassiola & Kovner(1993)) = b0 = E0(Eliasdottir (2007)) * (rs^2 - ra^2) / rs^2. So when s->\\infty and a->0, they are equivalent. In this implementation: diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py index 99d637963..93c26428a 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py @@ -5,7 +5,7 @@ from autogalaxy.profiles.mass.abstract.abstract import MassProfile -class dPIEP(MassProfile): +class dPIEPotential(MassProfile): def __init__( self, @@ -16,12 +16,12 @@ def __init__( b0: float = 1.0, ): """ - The dual Pseudo Isothermal Elliptical Potential (dPIEP) with pseudo-ellipticity on potential, based on the + The dual Pseudo Isothermal Elliptical Potential (dPIEPotential) with pseudo-ellipticity on potential, 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 dPIEP (Elliptical), this version + `kappa_scale`. Although originally called the dPIEPotential (Elliptical), this version lacks ellipticity, so the "E" may be a misnomer. The projected surface mass density is given by: @@ -69,7 +69,7 @@ def _ellip(self): def _deflection_angle(self, radii): """ - For a circularly symmetric dPIEP profile, computes the magnitude of the deflection at each radius. + For a circularly symmetric dPIEPotential profile, computes the magnitude of the deflection at each radius. """ a, s = self.ra, self.rs radii = jnp.maximum(radii, 1e-8) @@ -166,7 +166,7 @@ def potential_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): return jnp.zeros(shape=grid.shape[0]) -class dPIEPSph(dPIEP): +class dPIEPotentialSph(dPIEPotential): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -175,12 +175,12 @@ def __init__( b0: float = 1.0, ): """ - The dual Pseudo-Isothermal mass profile (dPIEP) without ellipticity, based on the + The dual Pseudo-Isothermal mass profile (dPIEPotential) 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 dPIEP (Elliptical), this version + `kappa_scale`. Although originally called the dPIEPotential (Elliptical), this version lacks ellipticity, so the "E" may be a misnomer. The projected surface mass density is given by: @@ -199,7 +199,7 @@ def __init__( `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. - Note: This dPIEPSph should be the same with dPIEMDSph for their same mathamatical formulations. + Note: This dPIEPotentialSph should be the same with dPIEMDSph for their same mathamatical formulations. Parameters ---------- diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py index dfe23df84..486ac1129 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py @@ -6,38 +6,38 @@ def test__deflections_yx_2d_from(): - mp = ag.mp.dPIEPSph(centre=(-0.7, 0.5), b0=5.2, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotentialSph(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) - mp = ag.mp.dPIEPSph(centre=(-0.1, 0.1), b0=20.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotentialSph(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) - mp = ag.mp.dPIEP(centre=(0, 0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential(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.186341843, 1e-3) assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) - mp = ag.mp.dPIEP(centre=(0, 0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential(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.186341843, 1e-3) assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) - elliptical = ag.mp.dPIEP( + elliptical = ag.mp.dPIEPotential( centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.0, ra=2.0, rs=3.0 ) - spherical = ag.mp.dPIEPSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) + spherical = ag.mp.dPIEPotentialSph(centre=(1.1, 1.1), b0=12.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 @@ -48,25 +48,25 @@ def test__convergence_2d_from(): # eta = 1.0 # kappa = 0.5 * 1.0 ** 1.0 - mp = ag.mp.dPIEPSph(centre=(0.0, 0.0), b0=8.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotentialSph(centre=(0.0, 0.0), b0=8.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.dPIEP(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=4.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=4.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.dPIEP(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=8.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=8.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.dPIEP( + mp = ag.mp.dPIEPotential( centre=(0.0, 0.0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0 ) @@ -74,10 +74,10 @@ def test__convergence_2d_from(): assert convergence == pytest.approx(0.87182837, 1e-3) - elliptical = ag.mp.dPIEP( + elliptical = ag.mp.dPIEPotential( centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.0, ra=2.0, rs=3.0 ) - spherical = ag.mp.dPIEPSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) + spherical = ag.mp.dPIEPotentialSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) assert elliptical.convergence_2d_from(grid=grid).array == pytest.approx( spherical.convergence_2d_from(grid=grid).array, 1e-4 From 132baebe436683ec82c87c1b5c20653594b019ce Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:15:49 +0100 Subject: [PATCH 5/9] dPIEMD -> dPIEMass --- .../total/dual_psuedo_isothermal_mass.yaml | 126 ++++++++++++++++++ autogalaxy/profiles/mass/__init__.py | 6 +- autogalaxy/profiles/mass/total/__init__.py | 2 +- .../mass/total/dual_pseudo_isothermal_mass.py | 36 ++--- .../total/dual_pseudo_isothermal_potential.py | 2 +- 5 files changed, 149 insertions(+), 23 deletions(-) create mode 100644 autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml diff --git a/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml b/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml new file mode 100644 index 000000000..7b5298aca --- /dev/null +++ b/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml @@ -0,0 +1,126 @@ +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 diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index 5aaa45e70..a4b358d0d 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -1,9 +1,9 @@ from .abstract.abstract import MassProfile from .point import PointMass, SMBH, SMBHBinary from .total import ( - dPIEMD, - dPIEMDSph, - PIEMD, + dPIEMass, + dPIEMassSph, + PIEMass, dPIEPotential, dPIEPotentialSph, PowerLawCore, diff --git a/autogalaxy/profiles/mass/total/__init__.py b/autogalaxy/profiles/mass/total/__init__.py index 1b0ecc88b..88ad4d76f 100644 --- a/autogalaxy/profiles/mass/total/__init__.py +++ b/autogalaxy/profiles/mass/total/__init__.py @@ -1,5 +1,5 @@ from .dual_pseudo_isothermal_potential import dPIEPotential, dPIEPotentialSph -from .dual_pseudo_isothermal_mass import PIEMD, dPIEMD, dPIEMDSph +from .dual_pseudo_isothermal_mass import PIEMass, dPIEMass, dPIEMassSph 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 be2ccd09f..d697aed35 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -4,13 +4,13 @@ import autoarray as aa from autogalaxy.profiles.mass.abstract.abstract import MassProfile -# Within this profile family, PIEMD, dPIEMD, and dPIEMDSph are directly ported from Lenstool's C code, and have been thoroughly annotated and adapted for PyAutoLens. +# 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. def _ci05(x, y, eps, rcore): """ - Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / E0 for PIEMD at any positions (x,y), + Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / E0 for PIEMass at any positions (x,y), see Kassiola & Kovner(1993) Eq. 4.1.2, which is the integral of Eq. 2.3.8. Note here b0(or called E0) is out of the `_ci05`. @@ -62,12 +62,12 @@ def _ci05(x, y, eps, rcore): def _ci05f(x, y, eps, rcore, rcut): """ - Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / (b0 * ra / (rs - ra)) for dPIEMD at any positions (x,y), + Returns the first derivatives of the lens potential as complex number I'* = (∂ψ/∂x + i ∂ψ/∂y) / (b0 * ra / (rs - ra)) for dPIEMass at any positions (x,y), which is the integral of Eq. 2.3.8 in Kassiola & Kovner(1993). - Note here (b0 * ra / (rs - ra)) is out of the `_ci05f`. The only difference of integral of Eq. 2.3.8 between dPIEMD and PIEMD is the \\kappa: - \\kappa(r_{em})_{dPIEMD} = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}). - I*_{dPIEMD} = ra / (rs - ra) * (I*_{PIEMD}(ra) - I*_{PIEMD}(ra)) + Note here (b0 * ra / (rs - ra)) is out of the `_ci05f`. The only difference of integral of Eq. 2.3.8 between dPIEMass and PIEMass is the \\kappa: + \\kappa(r_{em})_{dPIEMass} = rs / (rs - ra) * (\\kappa_{PIEMass,ra} - \\kappa_{PIEMass,rs}). + I*_{dPIEMass} = ra / (rs - ra) * (I*_{PIEMass}(ra) - I*_{PIEMass}(ra)) Parameters ---------- @@ -135,7 +135,7 @@ def _ci05f(x, y, eps, rcore, rcut): def _mdci05(x, y, eps, rcore, b0): """ - Returns the second derivatives (Hessian matrix) of the lens potential as complex number for PIEMD at any positions (x,y): + Returns the second derivatives (Hessian matrix) of the lens potential as complex number for PIEMass at any positions (x,y): ∂²ψ/∂x² = Re(∂I*/∂x), ∂²ψ/∂y² = Im(∂I*/∂y), ∂²ψ/∂x∂y = ∂²ψ/∂y∂x = Im(∂I*/∂x) = Re(∂I*/∂y) see Kassiola & Kovner(1993) Eq. 4.1.4. @@ -221,7 +221,7 @@ def _mdci05(x, y, eps, rcore, b0): return a, b, c, d -class PIEMD(MassProfile): +class PIEMass(MassProfile): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -230,7 +230,7 @@ def __init__( b0: float = 0.1, ): """ - The Pseudo Isothermal Elliptical Mass Distribution(PIEMD) profiles, based on the formulaiton from + The Pseudo Isothermal Elliptical Mass Distribution(PIEMass) profiles, based on the formulaiton from Kassiola & Kovner(1993) https://articles.adsabs.harvard.edu/pdf/1993ApJ...417..450K. This profile is ported from Lenstool's C code, which has the same formulation. @@ -325,7 +325,7 @@ def analytical_magnification_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs) return aa.Array2D(values=1.0 / det_A, mask=grid.mask) -class dPIEMD(MassProfile): +class dPIEMass(MassProfile): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -335,7 +335,7 @@ def __init__( b0: float = 0.1, ): """ - The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles, which is a *two component PIEMD* with both a core radius and a truncation radius, + The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMass) profiles, which is a *two component PIEMass* with both a core radius and a truncation radius, see Eliasdottir (2007): https://arxiv.org/abs/0710.5636 This profile is ported from Lenstool's C code, which has the same formulation. @@ -343,8 +343,8 @@ def __init__( and \\rho \~ r^{-4} in the outer parts: \\rho \\propto [(ra^2 + R^2) (rs^2 + R^2)]^{-1} - The convergence is given by two PIEMD with core radius ra and rs: - \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) + The convergence is given by two PIEMass with core radius ra and rs: + \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMass,ra} - \\kappa_{PIEMass,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEPotential}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. @@ -412,7 +412,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): """ Calculate the hessian matrix on a grid of (y,x) arc-second coordinates. - Hessian_dPIEMD = rs * (rs - ra) * ( Hessian_PIEMD(ra) - Hessian_PIEMD(rs)) + Hessian_dPIEMass = rs * (rs - ra) * ( Hessian_PIEMass(ra) - Hessian_PIEMass(rs)) Parameters ---------- @@ -451,7 +451,7 @@ def analytical_magnification_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs) return aa.Array2D(values=1.0 / det_A, mask=grid.mask) -class dPIEMDSph(dPIEMD): +class dPIEMassSph(dPIEMass): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), @@ -460,7 +460,7 @@ def __init__( b0: float = 1.0, ): """ - The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMD) profiles without ellipticity, which is a *two component PIEMD* with both a core radius and a truncation radius, + The dual Pseudo Isothermal Elliptical Mass Distribution(dPIEMass) profiles without ellipticity, which is a *two component PIEMass* with both a core radius and a truncation radius, see Eliasdottir (2007): https://arxiv.org/abs/0710.5636 This profile is ported from Lenstool's C code, which has the same formulation. @@ -468,8 +468,8 @@ def __init__( and \\rho \~ r^{-4} in the outer parts: \\rho \\propto [(ra^2 + R^2) (rs^2 + R^2)]^{-1} - The convergence is given by two PIEMD with core radius ra and rs: - \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMD,ra} - \\kappa_{PIEMD,rs}) + The convergence is given by two PIEMass with core radius ra and rs: + \\kappa(r_{em}) = rs / (rs - ra) * (\\kappa_{PIEMass,ra} - \\kappa_{PIEMass,rs}) = b_0 / 2 * rs / (rs - ra) * ( \\frac{1}{\\sqrt{ ra^2 + r_{em}^2}} - \\frac{1}{\\sqrt{ rs^2 + r_{em}^2}} ) where r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2. Note in Eliasdottir (2007), E0 = 6\\pi * \\sigma_{dPIEPotential}^2 / c^2 * (D_{LS} / D_{S}). Eliasdottir's E0 is not the same as E0 in Kassiola & Kovner(1993) which is also b0. diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py index 93c26428a..798b36989 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py @@ -199,7 +199,7 @@ def __init__( `b0` is not in the Intermediate-Axis-Convention for its r_{em}^2 = x^2 / (1 + \\epsilon)^2 + y^2 / (1 - \\epsilon)^2 Credit: Jackson O'Donnell for implementing this profile in PyAutoLens. - Note: This dPIEPotentialSph should be the same with dPIEMDSph for their same mathamatical formulations. + Note: This dPIEPotentialSph should be the same with dPIEMassSph for their same mathamatical formulations. Parameters ---------- From 884872c5f726f8a23b5d528e24e8f6a343b9b803 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:21:40 +0100 Subject: [PATCH 6/9] unit tests added --- .../mass/total/dual_pseudo_isothermal_mass.py | 24 ++++++------ .../total/test_dual_pseudo_isothermal_mass.py | 39 +++++++++++++++++++ .../test_dual_pseudo_isothermal_potential.py | 7 ---- 3 files changed, 51 insertions(+), 19 deletions(-) create mode 100644 test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index d697aed35..044b458d9 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -281,7 +281,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): """ ellip = self._ellip() factor = self.b0 - zis = _ci05(x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra) + zis = _ci05(x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra) # This is in axes aligned to the major/minor axis deflection_x = zis.real @@ -309,7 +309,7 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): ellip = self._ellip() hessian_xx, hessian_xy, hessian_yx, hessian_yy = _mdci05( - x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, b0=self.b0 + x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, b0=self.b0 ) return hessian_yy, hessian_xy, hessian_yx, hessian_xx @@ -396,7 +396,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): """ ellip = self._ellip() factor = self.b0 * self.rs / (self.rs - self.ra) - zis = _ci05f(x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, rcut=self.rs) + zis = _ci05f(x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, rcut=self.rs) # This is in axes aligned to the major/minor axis deflection_x = zis.real @@ -426,10 +426,10 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): t05 = self.rs / (self.rs - self.ra) g05c_a, g05c_b, g05c_c, g05c_d = _mdci05( - x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.ra, b0=self.b0 + x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, b0=self.b0 ) g05cut_a, g05cut_b, g05cut_c, g05cut_d = _mdci05( - x=grid[:, 1], y=grid[:, 0], eps=ellip, rcore=self.rs, b0=self.b0 + x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.rs, b0=self.b0 ) # Compute Hessian matrix components @@ -527,13 +527,13 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): s = self.rs # radii = self.radial_grid_from(grid=grid, **kwargs) # R2 = radii * radii - R2 = grid[:, 1] * grid[:, 1] + grid[:, 0] * grid[:, 0] + R2 = grid.array[:, 1] * grid.array[:, 1] + grid.array[:, 0] * grid.array[:, 0] factor = np.sqrt(R2 + a * a) - a - np.sqrt(R2 + s * s) + s factor *= self.b0 * s / (s - a) / R2 # This is in axes aligned to the major/minor axis - deflection_x = grid[:, 1] * factor - deflection_y = grid[:, 0] * factor + deflection_x = grid.array[:, 1] * factor + deflection_y = grid.array[:, 0] * factor # And here we convert back to the real axes return self.rotated_grid_from_reference_frame_from( @@ -580,14 +580,14 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): # ∂²ψ/∂²R = ∂(z/R)/∂R = (∂z/∂R * R - z * 1) / R^2 # = {( R^2 / √(R^2 + a^2)) - ( R^2 / √(R^2 + s^2)) - z} / R^2 # = p - R2 = grid[:, 1] * grid[:, 1] + grid[:, 0] * grid[:, 0] + R2 = grid.array[:, 1] * grid.array[:, 1] + grid.array[:, 0] * grid.array[:, 0] z = np.sqrt(R2 + a * a) - a - np.sqrt(R2 + s * s) + s p = (1.0 - a / np.sqrt(a * a + R2)) * a / R2 - ( 1.0 - s / np.sqrt(s * s + R2) ) * s / R2 - X = grid[:, 1] * grid[:, 1] / R2 # x^2 / R^2 - Y = grid[:, 0] * grid[:, 0] / R2 # y^2 / R^2 - XY = grid[:, 1] * grid[:, 0] / R2 # x*y / R^2 + X = grid.array[:, 1] * grid.array[:, 1] / R2 # x^2 / R^2 + Y = grid.array[:, 0] * grid.array[:, 0] / R2 # y^2 / R^2 + XY = grid.array[:, 1] * grid.array[:, 0] / R2 # x*y / R^2 # ∂²ψ/∂x² = ∂²ψ/∂R² * (∂R/∂x)² + ∂²R/∂x² * ∂ψ/∂R # = p * (x / R)^2 + y^2 / R^3 * z / R 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 new file mode 100644 index 000000000..ba1087083 --- /dev/null +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py @@ -0,0 +1,39 @@ +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.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) + + 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) + + # 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) + + elliptical = ag.mp.dPIEMass( + centre=(1.1, 1.1), ell_comps=(0.0, 0.0), 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-4 + ) diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py index 486ac1129..2b8d4174e 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py @@ -27,13 +27,6 @@ def test__deflections_yx_2d_from(): assert deflections[0, 0] == pytest.approx(0.186341843, 1e-3) assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) - mp = ag.mp.dPIEPotential(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.186341843, 1e-3) - assert deflections[0, 1] == pytest.approx(0.13176363087, 1e-3) - elliptical = ag.mp.dPIEPotential( centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.0, ra=2.0, rs=3.0 ) From bfdfccf909a7cd22ee14942202304cbd63a1116a Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 10:50:24 +0100 Subject: [PATCH 7/9] remove if statement to fix JAX --- ...al_mass.yaml => dual_pseudo_isothermal_mass.yaml} | 0 ...al.yaml => dual_pseudo_isothermal_potential.yaml} | 4 ++-- .../mass/total/dual_pseudo_isothermal_mass.py | 8 +++++++- .../mass/total/dual_pseudo_isothermal_potential.py | 9 +-------- .../mass/total/test_dual_pseudo_isothermal_mass.py | 4 +++- .../total/test_dual_pseudo_isothermal_potential.py | 12 +++++++++--- 6 files changed, 22 insertions(+), 15 deletions(-) rename autogalaxy/config/priors/mass/total/{dual_psuedo_isothermal_mass.yaml => dual_pseudo_isothermal_mass.yaml} (100%) rename autogalaxy/config/priors/mass/total/{dual_psuedo_isothermal_potential.yaml => dual_pseudo_isothermal_potential.yaml} (93%) diff --git a/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml b/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml similarity index 100% rename from autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_mass.yaml rename to autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml diff --git a/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml b/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_potential.yaml similarity index 93% rename from autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml rename to autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_potential.yaml index 4fe6d9e84..621ff3607 100644 --- a/autogalaxy/config/priors/mass/total/dual_psuedo_isothermal_potential.yaml +++ b/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_potential.yaml @@ -1,4 +1,4 @@ -dPIEP: +dPIEPotential: centre_0: type: Gaussian mean: 0.0 @@ -73,7 +73,7 @@ dPIEP: limits: lower: 0.0 upper: inf -dPIEPSph: +dPIEPotentialSph: centre_0: type: Gaussian mean: 0.0 diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index 044b458d9..ae64d69b9 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -396,7 +396,13 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): """ ellip = self._ellip() factor = self.b0 * self.rs / (self.rs - self.ra) - zis = _ci05f(x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, rcut=self.rs) + zis = _ci05f( + x=grid.array[:, 1], + y=grid.array[:, 0], + eps=ellip, + rcore=self.ra, + rcut=self.rs, + ) # This is in axes aligned to the major/minor axis deflection_x = zis.real diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py index 798b36989..d8d02ce52 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_potential.py @@ -55,9 +55,6 @@ def __init__( """ super().__init__(centre=centre, ell_comps=ell_comps) - if ra > rs: - ra, rs = rs, ra - self.ra = ra self.rs = rs self.b0 = b0 @@ -65,7 +62,7 @@ def __init__( def _ellip(self): ellip = jnp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) MAX_ELLIP = 0.99999 - return min(ellip, MAX_ELLIP) + return jnp.min(jnp.array([ellip, MAX_ELLIP])) def _deflection_angle(self, radii): """ @@ -213,10 +210,6 @@ def __init__( The lens strength in arcseconds. """ - # Ensure rs > ra (things will probably break otherwise) - if ra > rs: - ra, rs = rs, ra - super().__init__(centre=centre, ell_comps=(0.0, 0.0)) self.ra = ra 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 ba1087083..9be3c6e23 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 @@ -22,7 +22,9 @@ def test__deflections_yx_2d_from(): # 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) + 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]])) diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py index 2b8d4174e..ea43e7c06 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_potential.py @@ -20,7 +20,9 @@ def test__deflections_yx_2d_from(): assert deflections[0, 0] == pytest.approx(1.4212977207, 1e-4) assert deflections[0, 1] == pytest.approx(0.308977765378, 1e-4) - mp = ag.mp.dPIEPotential(centre=(0, 0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential( + 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]])) @@ -47,13 +49,17 @@ def test__convergence_2d_from(): assert convergence == pytest.approx(1.57182995, 1e-3) - mp = ag.mp.dPIEPotential(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=4.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential( + centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=4.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.dPIEPotential(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=8.0, ra=2.0, rs=3.0) + mp = ag.mp.dPIEPotential( + centre=(0.0, 0.0), ell_comps=(0.0, 0.0), b0=8.0, ra=2.0, rs=3.0 + ) convergence = mp.convergence_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0]])) From 9e72e011cd69e19202bc93714548e9f31e75ec39 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 11:02:04 +0100 Subject: [PATCH 8/9] update unitt est --- .../mass/total/dual_pseudo_isothermal_mass.py | 89 +++++++++---------- .../total/test_dual_pseudo_isothermal_mass.py | 4 +- 2 files changed, 43 insertions(+), 50 deletions(-) diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index ae64d69b9..ca1fc1519 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -1,4 +1,5 @@ from typing import Tuple +import jax.numpy as jnp import numpy as np import autoarray as aa @@ -25,9 +26,7 @@ def _ci05(x, y, eps, rcore): complex The value of the I'* term. """ - if eps < 1e-10: - eps = 1e-10 - sqe = np.sqrt(eps) + sqe = jnp.sqrt(eps) axis_ratio = (1.0 - eps) / (1.0 + eps) cxro = (1.0 + eps) * (1.0 + eps) cyro = (1.0 - eps) * (1.0 - eps) @@ -35,24 +34,24 @@ def _ci05(x, y, eps, rcore): ##### I'* = zres = zci * ln(zis) = zci * ln(znum / zden), see Eq. 4.1.2 ##### # Define intermediate complex variables - zci = np.complex128(complex(0.0, -0.5 * (1.0 - eps * eps) / sqe)) - znum = np.complex128( + zci = jnp.array(0.0 + 1j * (-0.5 * (1.0 - eps * eps) / sqe), dtype=jnp.complex128) + znum = jnp.complex128( axis_ratio * x - + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio) + + 1j * (2.0 * sqe * jnp.sqrt(rcore * rcore + rem2) - y / axis_ratio) ) - zden = np.complex128(x + 1j * (2.0 * rcore * sqe - y)) + zden = jnp.complex128(x + 1j * (2.0 * rcore * sqe - y)) # zis = znum / zden = (a+bi)/(c+di) = [(ac+bd)+(bc-ad i)] / (c^2+d^2) norm = zden.real * zden.real + zden.imag * zden.imag # |zden|^2 zis_re = (znum.real * zden.real + znum.imag * zden.imag) / norm zis_im = (znum.imag * zden.real - znum.real * zden.imag) / norm - zis = np.complex128(zis_re + 1j * zis_im) + zis = jnp.complex128(zis_re + 1j * zis_im) # ln(zis) = ln(|zis|) + i*Arg(zis) - zis_mag = np.abs(zis) - zis_re = np.log(zis_mag) - zis_im = np.angle(zis) - zis = np.complex128(zis_re + 1j * zis_im) + zis_mag = jnp.abs(zis) + zis_re = jnp.log(zis_mag) + zis_im = jnp.angle(zis) + zis = jnp.complex128(zis_re + 1j * zis_im) # I'* = zres = zci * ln(zis) zres = zci * zis @@ -82,9 +81,7 @@ def _ci05f(x, y, eps, rcore, rcut): complex The value of the I'* term. """ - if eps < 1e-10: - eps = 1e-10 - sqe = np.sqrt(eps) + sqe = jnp.sqrt(eps) axis_ratio = (1.0 - eps) / (1.0 + eps) cxro = (1.0 + eps) * (1.0 + eps) cyro = (1.0 - eps) * (1.0 - eps) @@ -93,16 +90,17 @@ def _ci05f(x, y, eps, rcore, rcut): ##### I'* = zres_rc - zres_rcut = zci * ln(zis_rc / zis_rcut) = zci * ln((znum_rc / zden_rc) / (znum_rcut / zden_rcut)) ##### # Define intermediate complex variables - zci = np.complex128(complex(0.0, -0.5 * (1.0 - eps * eps) / sqe)) - znum_rc = np.complex128( + zci = jnp.array(0.0 + 1j * (-0.5 * (1.0 - eps * eps) / sqe), dtype=jnp.complex128) + znum_rc = jnp.complex128( axis_ratio * x - + 1j * (2.0 * sqe * np.sqrt(rcore * rcore + rem2) - y / axis_ratio) + + 1j * (2.0 * sqe * jnp.sqrt(rcore * rcore + rem2) - y / axis_ratio) ) # a + bi - zden_rc = np.complex128(x + 1j * (2.0 * rcore * sqe - y)) # c + di - znum_rcut = np.complex128( - axis_ratio * x + 1j * (2.0 * sqe * np.sqrt(rcut * rcut + rem2) - y / axis_ratio) + zden_rc = jnp.complex128(x + 1j * (2.0 * rcore * sqe - y)) # c + di + znum_rcut = jnp.complex128( + axis_ratio * x + + 1j * (2.0 * sqe * jnp.sqrt(rcut * rcut + rem2) - y / axis_ratio) ) # a + ei - zden_rcut = np.complex128(x + 1j * (2.0 * rcut * sqe - y)) # c + fi + zden_rcut = jnp.complex128(x + 1j * (2.0 * rcut * sqe - y)) # c + fi # zis_rc = znum_rc / zden_rc = (a+bi)/(c+di) # zis_rcut = znum_rcut / zden_rcut = (a+ei)/(c+fi) @@ -119,13 +117,13 @@ def _ci05f(x, y, eps, rcore, rcut): norm = cc * cc + dd * dd aaa = (aa * cc + bb * dd) / norm bbb = (bb * cc - aa * dd) / norm - zis_tot = np.complex128(aaa + 1j * bbb) + zis_tot = jnp.complex128(aaa + 1j * bbb) # ln(zis_tot) = ln(|zis_tot|) + i*Arg(zis_tot) - zis_tot_mag = np.abs(zis_tot) - zr_re = np.log(zis_tot_mag) - zr_im = np.angle(zis_tot) - zr = np.complex128(zr_re + 1j * zr_im) + zis_tot_mag = jnp.abs(zis_tot) + zr_re = jnp.log(zis_tot_mag) + zr_im = jnp.angle(zis_tot) + zr = jnp.complex128(zr_re + 1j * zr_im) # I'* = zci * ln(zis_tot) zres = zci * zr @@ -150,20 +148,18 @@ def _mdci05(x, y, eps, rcore, b0): complex The value of the I'* term. """ - if eps < 1e-10: - eps = 1e-10 # Calculate intermediate values # I*(x,y) = b0 * ci * (-i) * (ln{ q * x + (2.0 * sqe * wrem - y * 1/q )*i} - ln{ x + (2.0 * rcore * sqe - y)*i}) # = b0 * ci * (-i) * (ln{ q * x + num1*i} - ln{ x + num2*i}) # = b0 * ci * (-i) * (ln{u(x,y)} - ln{v(x,y)}) - sqe = np.sqrt(eps) + sqe = jnp.sqrt(eps) axis_ratio = (1.0 - eps) / (1.0 + eps) axis_ratio_inv = 1.0 / axis_ratio cxro = (1.0 + eps) * (1.0 + eps) cyro = (1.0 - eps) * (1.0 - eps) ci = 0.5 * (1.0 - eps * eps) / sqe - wrem = np.sqrt(rcore * rcore + x * x / cxro + y * y / cyro) # √(w(x,y)) + wrem = jnp.sqrt(rcore * rcore + x * x / cxro + y * y / cyro) # √(w(x,y)) num1 = 2.0 * sqe * wrem - y * axis_ratio_inv den1 = axis_ratio * axis_ratio * x * x + num1 * num1 # |q * x + num1*i|^2 num2 = 2.0 * rcore * sqe - y @@ -264,9 +260,9 @@ def __init__( self.b0 = b0 def _ellip(self): - ellip = np.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) + ellip = jnp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) MAX_ELLIP = 0.99999 - return min(ellip, MAX_ELLIP) + return jnp.min(jnp.array([ellip, MAX_ELLIP])) @aa.grid_dec.to_vector_yx @aa.grid_dec.transform @@ -289,7 +285,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # And here we convert back to the real axes return self.rotated_grid_from_reference_frame_from( - grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), + grid=jnp.multiply(factor, jnp.vstack((deflection_y, deflection_x)).T), **kwargs, ) @@ -303,7 +299,7 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ - grid = np.asarray(grid) + grid = jnp.asarray(grid) if grid.ndim != 2 or grid.shape[1] != 2: raise ValueError("Grid must be a 2D array with shape (n, 2)") ellip = self._ellip() @@ -371,17 +367,14 @@ def __init__( """ super().__init__(centre=centre, ell_comps=ell_comps) - if ra > rs: - ra, rs = rs, ra - self.ra = ra self.rs = rs self.b0 = b0 def _ellip(self): - ellip = np.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) + ellip = jnp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) MAX_ELLIP = 0.99999 - return min(ellip, MAX_ELLIP) + return jnp.min(jnp.array([ellip, MAX_ELLIP])) @aa.grid_dec.to_vector_yx @aa.grid_dec.transform @@ -410,7 +403,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # And here we convert back to the real axes return self.rotated_grid_from_reference_frame_from( - grid=np.multiply(factor, np.vstack((deflection_y, deflection_x)).T), + grid=jnp.multiply(factor, jnp.vstack((deflection_y, deflection_x)).T), **kwargs, ) @@ -425,7 +418,7 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ - grid = np.asarray(grid) + grid = jnp.asarray(grid) if grid.ndim != 2 or grid.shape[1] != 2: raise ValueError("Grid must be a 2D array with shape (n, 2)") ellip = self._ellip() @@ -534,7 +527,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # radii = self.radial_grid_from(grid=grid, **kwargs) # R2 = radii * radii R2 = grid.array[:, 1] * grid.array[:, 1] + grid.array[:, 0] * grid.array[:, 0] - factor = np.sqrt(R2 + a * a) - a - np.sqrt(R2 + s * s) + s + factor = jnp.sqrt(R2 + a * a) - a - jnp.sqrt(R2 + s * s) + s factor *= self.b0 * s / (s - a) / R2 # This is in axes aligned to the major/minor axis @@ -543,7 +536,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs): # 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 + grid=jnp.multiply(1.0, jnp.vstack((deflection_y, deflection_x)).T), **kwargs ) @aa.grid_dec.transform @@ -561,7 +554,7 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ - grid = np.asarray(grid) + grid = jnp.asarray(grid) if grid.ndim != 2 or grid.shape[1] != 2: raise ValueError("Grid must be a 2D array with shape (n, 2)") @@ -587,9 +580,9 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", **kwargs): # = {( R^2 / √(R^2 + a^2)) - ( R^2 / √(R^2 + s^2)) - z} / R^2 # = p R2 = grid.array[:, 1] * grid.array[:, 1] + grid.array[:, 0] * grid.array[:, 0] - z = np.sqrt(R2 + a * a) - a - np.sqrt(R2 + s * s) + s - p = (1.0 - a / np.sqrt(a * a + R2)) * a / R2 - ( - 1.0 - s / np.sqrt(s * s + R2) + z = jnp.sqrt(R2 + a * a) - a - jnp.sqrt(R2 + s * s) + s + p = (1.0 - a / jnp.sqrt(a * a + R2)) * a / R2 - ( + 1.0 - s / jnp.sqrt(s * s + R2) ) * s / R2 X = grid.array[:, 1] * grid.array[:, 1] / R2 # x^2 / R^2 Y = grid.array[:, 0] * grid.array[:, 0] / R2 # y^2 / R^2 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 9be3c6e23..ef5b441c0 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 @@ -32,10 +32,10 @@ def test__deflections_yx_2d_from(): assert deflections[0, 1] == pytest.approx(0.10753914, 1e-3) elliptical = ag.mp.dPIEMass( - centre=(1.1, 1.1), ell_comps=(0.0, 0.0), b0=12.0, ra=2.0, rs=3.0 + 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-4 + spherical.deflections_yx_2d_from(grid=grid).array, 1e-1 ) From e9c369af6145f12e8c4ae27a4f76736019c9a0d8 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 15 Oct 2025 11:06:55 +0100 Subject: [PATCH 9/9] JAX conversion complete and works --- autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index ca1fc1519..0211aa8ce 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -494,8 +494,7 @@ def __init__( The lens strength in arcseconds. """ super().__init__(centre=centre, ell_comps=(0.0, 0.0)) - if ra > rs: - ra, rs = rs, ra + self.ra = ra self.rs = rs self.b0 = b0