From f3cdf6ff43c898f50095398b1e7f1cf80519f59b Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 29 Oct 2026 21:25:59 +0000 Subject: [PATCH] fix: lmp.Gaussian potential as the line integral of its analytic deflection The potential decomposed a single Gaussian into 100 Gaussians and was wrong: its gradient did not reproduce the profile's deflections (sign errors for elliptical Gaussians). It is now psi = int_0^1 alpha(t r).r dt over the 20-node Gauss-Legendre rule, using the same in-frame analytic deflection as deflections_2d_via_analytic_from (one Faddeeva call per node). Gradient matches the deflections to 1e-8 and lap(psi) = 2 kappa; a 40-member MGE basis potential on 100x100 drops from >10 min to ~2 s. Closes #680 Co-Authored-By: Claude Opus 5.5 --- autogalaxy/profiles/mass/stellar/gaussian.py | 44 +++++++++++++--- .../profiles/mass/stellar/test_gaussian.py | 52 +++++++++++++++++++ 2 files changed, 88 insertions(+), 8 deletions(-) diff --git a/autogalaxy/profiles/mass/stellar/gaussian.py b/autogalaxy/profiles/mass/stellar/gaussian.py index e2b58bd3..eb458c08 100644 --- a/autogalaxy/profiles/mass/stellar/gaussian.py +++ b/autogalaxy/profiles/mass/stellar/gaussian.py @@ -94,6 +94,16 @@ def deflections_2d_via_analytic_from( """ + return self._deflections_in_frame_from(grid=grid, xp=xp) + + def _deflections_in_frame_from(self, grid: aa.type.Grid2DLike, xp=np): + """ + The (y,x) deflection angles on a ``grid`` already in the profile's reference + frame (centred and rotated), without the transform decorators. + + Shared by ``deflections_2d_via_analytic_from`` and ``potential_2d_from`` so the + potential is the exact line integral of the same deflection field. + """ if _is_circular(self.ell_comps): # A circular Gaussian has an exact real radial form -- take it instead of # the complex Faddeeva evaluated at the q = 0.9999 clamp. The predicate @@ -141,14 +151,32 @@ def convergence_func(self, grid_radius: float, xp=np) -> float: @aa.decorators.to_array @aa.decorators.transform def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): - radii_min = self.sigma / 100.0 - radii_max = self.sigma * 20.0 - sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 100)) - mge_decomp = MGEDecomposer(mass_profile=self) - return mge_decomp.potential_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="circularised", three_D=False, - ) + """ + The lensing potential, as the line integral of the analytic deflection field + from the centre to each grid point, + + psi(x, y) = integral_0^1 alpha(t x, t y) . (x, y) dt, + + with the 20-node Gauss-Legendre rule of ``MGEDecomposer``. A Gaussian is its own + one-term MGE, so this needs one Faddeeva evaluation per node rather than a + decomposition into many Gaussians. The potential is zero at the centre. + """ + y = xp.asarray(grid.array[:, 0], dtype=xp.float64) + x = xp.asarray(grid.array[:, 1], dtype=xp.float64) + + potential = xp.zeros(y.shape[0], dtype=xp.float64) + + for node, weight in zip( + MGEDecomposer._GL_NODES_20, MGEDecomposer._GL_WEIGHTS_20 + ): + t = 0.5 * (1.0 + float(node)) + scaled_grid = aa.Grid2DIrregular(values=xp.stack([t * y, t * x], axis=-1)) + deflections = self._deflections_in_frame_from(grid=scaled_grid, xp=xp) + potential = potential + 0.5 * float(weight) * ( + deflections[:, 0] * y + deflections[:, 1] * x + ) + + return potential def image_2d_via_radii_from(self, grid_radii: np.ndarray, xp=np): """Calculate the intensity of the Gaussian light profile on a grid of radial coordinates. diff --git a/test_autogalaxy/profiles/mass/stellar/test_gaussian.py b/test_autogalaxy/profiles/mass/stellar/test_gaussian.py index 8cc3ea6c..32cc8eef 100644 --- a/test_autogalaxy/profiles/mass/stellar/test_gaussian.py +++ b/test_autogalaxy/profiles/mass/stellar/test_gaussian.py @@ -347,3 +347,55 @@ def test__deflections_yx_2d_from__spherical_case__is_radial_and_matches_elliptic deflections_ell = np.asarray(gaussian_ell.deflections_yx_2d_from(grid=grid).array) assert deflections_sph == pytest.approx(deflections_ell, rel=1.0e-5) + + +@pytest.mark.parametrize("ell_comps", [(0.0, 0.0), (0.2, -0.1), (0.43, 0.18)]) +def test__potential_2d_from__gradient_is_deflections(ell_comps): + gaussian = ag.mp.Gaussian( + centre=(0.05, -0.03), + ell_comps=ell_comps, + intensity=2.0, + sigma=0.4, + mass_to_light_ratio=1.5, + ) + points = np.array([[0.7, 0.3], [-0.5, 0.9], [1.2, -1.0], [0.2, 0.1]]) + step = 1e-4 + + def potential(values): + return np.asarray( + gaussian.potential_2d_from(grid=ag.Grid2DIrregular(values)).array + ) + + gradient_y = (potential(points + [step, 0]) - potential(points - [step, 0])) / ( + 2 * step + ) + gradient_x = (potential(points + [0, step]) - potential(points - [0, step])) / ( + 2 * step + ) + deflections = np.asarray( + gaussian.deflections_yx_2d_from(grid=ag.Grid2DIrregular(points)).array + ) + + assert np.c_[gradient_y, gradient_x] == pytest.approx(deflections, abs=1e-6) + + +def test__potential_2d_from__zero_at_centre_and_basis_is_fast(): + gaussian = ag.mp.Gaussian( + centre=(0.1, 0.2), ell_comps=(0.3, 0.1), intensity=1.0, sigma=0.5 + ) + assert gaussian.potential_2d_from(grid=ag.Grid2DIrregular([(0.1, 0.2)])).array[ + 0 + ] == pytest.approx(0.0, abs=1e-12) + + import time + + basis = ag.lp_basis.Basis( + profile_list=[ + ag.lmp.Gaussian(ell_comps=(0.43, 0.18), intensity=1.0, sigma=sigma) + for sigma in np.logspace(-2, 0.5, 40) + ] + ) + grid = ag.Grid2D.uniform(shape_native=(50, 50), pixel_scales=0.1) + start = time.time() + basis.potential_2d_from(grid=grid) + assert time.time() - start < 30.0