From 0b69d8027a1b32e4f722ae552d59c3e6807b0635 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 29 May 2026 09:14:53 +0100 Subject: [PATCH] fix(mass): convergence_func on PowerLawBroken, PowerLawMultipole, cNFW family Follow-up to #466 (dPIE convergence_func). Adds the abstract convergence_func hook to the three mass-profile families that define convergence_2d_from but never overrode convergence_func, so radial mass integration (mass_integral -> mass_angular_within_circle_from -> Einstein radius) and MGE three_D=False potential decomposition no longer raise NotImplementedError. PowerLawBroken's MGE potential is incompatible with its piecewise convergence, so potential_2d_from now raises an explanatory NotImplementedError. Co-Authored-By: Claude Opus 4.8 (1M context) --- autogalaxy/profiles/mass/dark/cnfw.py | 84 +++++++++++++++---- .../profiles/mass/total/power_law_broken.py | 70 +++++++++++----- .../mass/total/power_law_multipole.py | 14 ++++ .../profiles/mass/dark/test_cnfw.py | 50 +++++++++++ .../mass/total/test_power_law_broken.py | 61 ++++++++++++++ .../mass/total/test_power_law_multipole.py | 28 +++++++ 6 files changed, 269 insertions(+), 38 deletions(-) diff --git a/autogalaxy/profiles/mass/dark/cnfw.py b/autogalaxy/profiles/mass/dark/cnfw.py index 18639b486..a6882ba15 100644 --- a/autogalaxy/profiles/mass/dark/cnfw.py +++ b/autogalaxy/profiles/mass/dark/cnfw.py @@ -100,12 +100,12 @@ class cNFW(AbstractgNFW): """ def __init__( - self, - centre: Tuple[float, float] = (0.0, 0.0), - ell_comps: Tuple[float, float] = (0.0, 0.0), - kappa_s: float = 0.05, - scale_radius: float = 1.0, - core_radius: float = 0.01, + self, + centre: Tuple[float, float] = (0.0, 0.0), + ell_comps: Tuple[float, float] = (0.0, 0.0), + kappa_s: float = 0.05, + scale_radius: float = 1.0, + core_radius: float = 0.01, ): r""" Parameters @@ -129,7 +129,6 @@ def __init__( self.scale_radius = scale_radius self.core_radius = core_radius - def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): return self.deflections_2d_via_mge_from(grid=grid, xp=xp, **kwargs) @@ -145,8 +144,8 @@ def deflections_2d_via_mge_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs) grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention='major', - three_D=True + ellipticity_convention="major", + three_D=True, ) return deflections_via_mge @@ -163,6 +162,45 @@ def density_3d_func(self, r, xp=np): * (r.array + self.scale_radius) ** (-2.0) ) + def convergence_func(self, grid_radius, xp=np): + """ + Radial projected convergence kappa(r), reusing the MGE-of-3D-density decomposition + (the same machinery `convergence_2d_from` uses with `three_D=True`) evaluated at the + radial coordinate `grid_radius`. + + cNFW has no closed-form radial convergence helper, so this delegates to the MGE + Gaussian sum. This hook is reached by `MGEDecomposer.decompose_convergence_via_mge` + (`three_D=False`, not used by cNFW) and by radial mass integration (`mass_integral` + -> `mass_angular_within_circle_from` -> Einstein radius). + + The result is the q-independent radial profile (like `NFW.convergence_func`): + ellipticity is re-introduced by the MGE machinery elsewhere, so no `sigmas_factor` + rescale is applied (`sigmas_factor=1.0`). Verified to match `convergence_2d_from` for + the spherical case and to be q-independent for the elliptical case. + """ + radii = ( + grid_radius.array + if hasattr(grid_radius, "array") + else xp.asarray(grid_radius) + ) + # Track scalar input so we return a scalar (matching other `convergence_func` + # implementations). `mass_integral` -> scipy.quad feeds scalar radii and a length-1 + # array return would warn (and eventually error) on the array->scalar conversion. + scalar_input = xp.ndim(radii) == 0 + radii = xp.atleast_1d(radii) + radii_min = self.scale_radius / 1000.0 + radii_max = self.scale_radius * 200.0 + sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 20)) + mge_decomp = MGEDecomposer(mass_profile=self) + convergence = mge_decomp._convergence_2d_via_mge_from( + grid_radii=radii, + xp=xp, + sigma_log_list=sigmas, + three_D=True, + sigmas_factor=1.0, + ) + return convergence[0] if scalar_input else convergence + @aa.over_sample @aa.decorators.to_array @aa.decorators.transform @@ -172,8 +210,11 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 20)) mge_decomp = MGEDecomposer(mass_profile=self) return mge_decomp.convergence_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=True, + grid=grid, + xp=xp, + sigma_log_list=sigmas, + ellipticity_convention="major", + three_D=True, ) @aa.over_sample @@ -185,8 +226,11 @@ def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 20)) mge_decomp = MGEDecomposer(mass_profile=self) return mge_decomp.potential_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=True, + grid=grid, + xp=xp, + sigma_log_list=sigmas, + ellipticity_convention="major", + three_D=True, ) @@ -308,8 +352,11 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 20)) mge_decomp = MGEDecomposer(mass_profile=self) return mge_decomp.convergence_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=True, + grid=grid, + xp=xp, + sigma_log_list=sigmas, + ellipticity_convention="major", + three_D=True, ) @aa.over_sample @@ -321,6 +368,9 @@ def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 20)) mge_decomp = MGEDecomposer(mass_profile=self) return mge_decomp.potential_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=True, + grid=grid, + xp=xp, + sigma_log_list=sigmas, + ellipticity_convention="major", + three_D=True, ) diff --git a/autogalaxy/profiles/mass/total/power_law_broken.py b/autogalaxy/profiles/mass/total/power_law_broken.py index f686f99a9..ab87741cf 100644 --- a/autogalaxy/profiles/mass/total/power_law_broken.py +++ b/autogalaxy/profiles/mass/total/power_law_broken.py @@ -90,40 +90,68 @@ def __init__( else: self.kB = (2 - self.inner_slope) / (2 * self.nu**2) - @aa.over_sample - @aa.decorators.to_array - @aa.decorators.transform - def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): - """ - Returns the dimensionless density kappa=Sigma/Sigma_c (eq. 1) + def _convergence(self, radii, xp=np): """ + Returns the dimensionless density kappa=Sigma/Sigma_c (eq. 1) as a function of the + (elliptical or circular) radial coordinate `radii`. - # Ell radius - radius = xp.hypot(grid.array[:, 1] * self.axis_ratio(xp), grid.array[:, 0]) + Shared by `convergence_2d_from` (which passes the elliptical radius) and the + `convergence_func` hook (which passes the radial coordinate directly). `radii` is a + plain array/scalar (callers unwrap `aa.ArrayIrregular` to `.array` first), so the + boolean masks `(radii <= break_radius)` return raw numpy bool arrays. + """ # Inside break radius - kappa_inner = self.kB * (self.break_radius / radius) ** self.inner_slope + kappa_inner = self.kB * (self.break_radius / radii) ** self.inner_slope # Outside break radius - kappa_outer = self.kB * (self.break_radius / radius) ** self.outer_slope + kappa_outer = self.kB * (self.break_radius / radii) ** self.outer_slope - return kappa_inner * (radius <= self.break_radius) + kappa_outer * ( - radius > self.break_radius + return kappa_inner * (radii <= self.break_radius) + kappa_outer * ( + radii > self.break_radius ) + def convergence_func(self, grid_radius, xp=np): + # Unwrap `aa.ArrayIrregular` -> plain array so `_convergence` returns a plain array + # (matching e.g. `Isothermal.convergence_func`). `mass_integral` -> scipy.quad cannot + # consume an `aa.ArrayIrregular` return. + radii = grid_radius.array if hasattr(grid_radius, "array") else grid_radius + return self._convergence(radii, xp=xp) + @aa.over_sample @aa.decorators.to_array @aa.decorators.transform + def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): + """ + Returns the dimensionless density kappa=Sigma/Sigma_c (eq. 1) + """ + + # Ell radius + radius = xp.hypot(grid.array[:, 1] * self.axis_ratio(xp), grid.array[:, 0]) + + return self._convergence(radius, xp=xp) + def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): - from autogalaxy.profiles.mass.abstract.mge import MGEDecomposer - - radii_min = self.break_radius / 100.0 - radii_max = self.einstein_radius * 20.0 - sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 30)) - mge_decomp = MGEDecomposer(mass_profile=self) - return mge_decomp.potential_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=False, + """ + The lensing potential is not available for the broken power law. + + It would be computed by decomposing the projected convergence into Gaussians via + `potential_2d_via_mge_from` (`three_D=False`), but that decomposition integrates the + convergence along a *complex* contour and therefore requires an analytic convergence + profile. The broken power law's convergence is piecewise — a slope discontinuity at + `break_radius` — so it is non-analytic and the MGE potential is numerically invalid + (wrong by many orders of magnitude). A correct potential would integrate this + profile's own analytic deflection field (eq. 18-19) along radial lines, which is not + yet implemented. + + `convergence_2d_from`, `deflections_yx_2d_from`, and `convergence_func` (and hence the + Einstein-radius / enclosed-mass integrals) are all available and correct. + """ + raise NotImplementedError( + "PowerLawBroken.potential_2d_from is not implemented: the MGE potential " + "decomposition requires an analytic convergence profile, but the broken power " + "law's convergence is piecewise (a kink at break_radius). Use deflections / " + "convergence instead, or integrate the analytic deflections for the potential." ) @aa.decorators.to_vector_yx diff --git a/autogalaxy/profiles/mass/total/power_law_multipole.py b/autogalaxy/profiles/mass/total/power_law_multipole.py index 942e24606..c012e06f8 100644 --- a/autogalaxy/profiles/mass/total/power_law_multipole.py +++ b/autogalaxy/profiles/mass/total/power_law_multipole.py @@ -289,6 +289,20 @@ def convergence_2d_from( * xp.cos(self.m * (angle - angle_m)) ) + def convergence_func(self, grid_radius, xp=np): + """ + The radial (azimuthally-averaged) convergence of a pure multipole perturbation is + identically zero — the `cos(m(phi - phi_m))` term integrates to zero over angle for + `m >= 1`, so the multipole adds no net azimuthally-symmetric mass. + + This hook is reached only by radial integration (`mass_integral` -> + `mass_angular_within_circle_from`); returning the zero monopole means a pure + multipole correctly encloses zero net mass. `.array` is unwrapped first so the + return is a plain array rather than an `aa.ArrayIrregular`. + """ + radii = grid_radius.array if hasattr(grid_radius, "array") else grid_radius + return xp.zeros_like(radii) + @aa.decorators.to_array def potential_2d_from( self, grid: aa.type.Grid2DLike, xp=np, **kwargs diff --git a/test_autogalaxy/profiles/mass/dark/test_cnfw.py b/test_autogalaxy/profiles/mass/dark/test_cnfw.py index 1dca64900..c2788c15b 100644 --- a/test_autogalaxy/profiles/mass/dark/test_cnfw.py +++ b/test_autogalaxy/profiles/mass/dark/test_cnfw.py @@ -16,3 +16,53 @@ def test__deflections_yx_2d_from(): deflection_r = np.sqrt(deflection_2d[0, 0] ** 2 + deflection_2d[0, 1] ** 2) assert deflection_r == pytest.approx(0.006034319441107217, 1.0e-8) + + +def test__convergence_func__matches_convergence_2d_from_and_is_q_independent(): + """Regression: the cNFW family overrides the abstract `convergence_func` (reached by + radial mass integration / Einstein radius) by reusing the MGE-of-3D-density + decomposition evaluated radially. It is q-independent (like `NFW.convergence_func`): + ellipticity is re-introduced by the MGE machinery elsewhere, so no `sigmas_factor` + rescale is applied.""" + + kappa_s, scale_radius, core_radius = 0.2, 2.0, 0.1 + radii = np.array([0.5, 1.0, 2.0, 3.0]) + + sph = ag.mp.cNFWSph( + centre=(0.0, 0.0), + kappa_s=kappa_s, + scale_radius=scale_radius, + core_radius=core_radius, + ) + + # `convergence_2d_from` on a `Grid2DIrregular` short-circuits the `@over_sample` + # decorator, so points at (0, r) evaluate the radial profile exactly. + grid = ag.Grid2DIrregular([[0.0, r] for r in radii]) + truth = np.asarray(sph.convergence_2d_from(grid=grid).array) + + actual = np.asarray(sph.convergence_func(ag.ArrayIrregular(radii))) + assert actual == pytest.approx(truth, 1e-8) + + # An elliptical cNFW returns the SAME radial (q-independent) convergence; this guards + # the `sigmas_factor=1.0` choice (a sqrt(q) factor would make the elliptical case wrong). + ell = ag.mp.cNFW( + centre=(0.0, 0.0), + ell_comps=(0.0, 0.4), + kappa_s=kappa_s, + scale_radius=scale_radius, + core_radius=core_radius, + ) + actual_ell = np.asarray(ell.convergence_func(ag.ArrayIrregular(radii))) + assert actual_ell == pytest.approx(truth, 1e-8) + + +def test__convergence_func__mass_integral_runs(): + """The Einstein-radius / enclosed-mass integral (`mass_angular_within_circle_from`) + routes through `mass_integral` -> `convergence_func` via scipy.quad; before the + override this raised NotImplementedError.""" + + cnfw = ag.mp.cNFWSph( + centre=(0.0, 0.0), kappa_s=0.2, scale_radius=2.0, core_radius=0.1 + ) + + assert cnfw.mass_angular_within_circle_from(radius=1.0) > 0.0 diff --git a/test_autogalaxy/profiles/mass/total/test_power_law_broken.py b/test_autogalaxy/profiles/mass/total/test_power_law_broken.py index 8139ea065..537fabfc5 100644 --- a/test_autogalaxy/profiles/mass/total/test_power_law_broken.py +++ b/test_autogalaxy/profiles/mass/total/test_power_law_broken.py @@ -230,3 +230,64 @@ def test__deflections_yx_2d_from__compare_to_power_law__slope_24(): power_law_yx_ratio = deflections[0, 0] / deflections[0, 1] assert broken_yx_ratio == pytest.approx(power_law_yx_ratio, 1.0e-4) + + +def test__convergence_func__matches_private_helper(): + """Regression: PowerLawBroken must override the abstract `convergence_func` so + `MGEDecomposer.decompose_convergence_via_mge` (reached by `potential_2d_from`, + which uses the `three_D=False` MGE path) doesn't fall through to the abstract + NotImplementedError. The shim delegates to the `_convergence` radial helper that + `convergence_2d_from` already uses.""" + + mp = ag.mp.PowerLawBroken( + centre=(0.0, 0.0), + ell_comps=(0.1, 0.05), + einstein_radius=1.0, + inner_slope=1.5, + outer_slope=2.5, + break_radius=0.4, + ) + + # Scalar radius, inside and outside the break radius. + assert mp.convergence_func(0.2) == pytest.approx(mp._convergence(0.2), 1e-12) + assert mp.convergence_func(1.5) == pytest.approx(mp._convergence(1.5), 1e-12) + + # 1-D array of radii spanning the break: shape preserved, element-wise equality. + radii = np.array([0.1, 0.4, 0.5, 1.0, 2.5]) + expected = mp._convergence(radii) + actual = mp.convergence_func(radii) + assert actual.shape == radii.shape + assert actual == pytest.approx(expected, 1e-12) + + # PowerLawBrokenSph inherits the override from PowerLawBroken. + sph = ag.mp.PowerLawBrokenSph( + centre=(0.0, 0.0), + einstein_radius=1.0, + inner_slope=1.5, + outer_slope=2.5, + break_radius=0.4, + ) + assert sph.convergence_func(0.2) == pytest.approx(sph._convergence(0.2), 1e-12) + + +def test__potential_2d_from__raises_not_implemented(): + """The broken power law's convergence is piecewise (a slope kink at break_radius), which + is incompatible with the MGE potential decomposition (it integrates the convergence along + a complex contour and requires an analytic profile, producing values wrong by orders of + magnitude otherwise). `potential_2d_from` therefore raises a clear NotImplementedError + rather than returning numerically-invalid values. Convergence/deflections/convergence_func + remain available.""" + + mp = ag.mp.PowerLawBroken( + centre=(0.0, 0.0), + ell_comps=(0.1, 0.05), + einstein_radius=1.0, + inner_slope=1.5, + outer_slope=2.5, + break_radius=0.4, + ) + + grid = ag.Grid2D.uniform(shape_native=(5, 5), pixel_scales=0.2) + + with pytest.raises(NotImplementedError): + mp.potential_2d_from(grid=grid) diff --git a/test_autogalaxy/profiles/mass/total/test_power_law_multipole.py b/test_autogalaxy/profiles/mass/total/test_power_law_multipole.py index 284d43721..b568df7e6 100644 --- a/test_autogalaxy/profiles/mass/total/test_power_law_multipole.py +++ b/test_autogalaxy/profiles/mass/total/test_power_law_multipole.py @@ -75,3 +75,31 @@ def test__potential_2d_from(): potential = mp.potential_2d_from(grid=ag.Grid2DIrregular([[1.0, 0.0]])) assert potential[0] == pytest.approx(0.0, 1e-3) + + +def test__convergence_func__returns_zero_monopole(): + """Regression: PowerLawMultipole overrides the abstract `convergence_func` to return + the zero monopole. The multipole convergence is angle-dependent (cos(m(phi - phi_m))) + so its azimuthal average is identically zero — a pure multipole encloses zero net + azimuthally-symmetric mass. `convergence_func` is reached only by radial mass + integration (`mass_integral` -> `mass_angular_within_circle_from`).""" + + import numpy as np + + mp = ag.mp.PowerLawMultipole( + m=4, + centre=(0.0, 0.0), + einstein_radius=1.0, + slope=2.0, + multipole_comps=(0.1, 0.2), + ) + + assert mp.convergence_func(1.5) == pytest.approx(0.0, 1e-12) + + radii = np.array([0.1, 0.5, 1.0, 2.5]) + actual = mp.convergence_func(radii) + assert actual.shape == radii.shape + assert actual == pytest.approx(np.zeros_like(radii), 1e-12) + + # A pure multipole encloses zero net mass. + assert mp.mass_angular_within_circle_from(radius=2.0) == pytest.approx(0.0, 1e-9)