fix: lmp.Gaussian potential from its analytic deflection (correct + ~100x faster)

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.

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