fix(mass): convergence_func on PowerLawBroken, PowerLawMultipole, cNFW family

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,
         )

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

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

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

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)

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)