Declare the ell_comps saturation constraint on EllProfile

autogalaxy/convert.py

@@ -1,6 +1,21 @@
 import numpy as np
 from typing import Tuple, Optional
 
+# The ellipticity magnitude at which the axis-ratio conversion saturates. Beyond
+# it every magnitude maps onto the same axis ratio, so the likelihood is exactly
+# flat in the radial direction and carries no gradient there — a gradient search
+# that walks past this cannot walk back. `EllProfile.__model_constraint__`
+# reports the distance beyond it, which is how PyAutoFit's multi-start searches
+# count the lanes trapped in that region.
+#
+# Deliberately distinct from the `magnitude_squared >= 1.0` validity threshold in
+# `profiles/validate.py`. They answer different questions — this is where the
+# *gradient* dies, that is where the *geometry* stops meaning anything — and the
+# annulus between them is reachable: at magnitude 0.9995 the radial derivative is
+# already exactly zero while `validate_ell_comps` still calls the point valid.
+# Stated together here rather than left as unrelated literals in separate files.
+ELL_COMPS_MAGNITUDE_CLAMP = 0.999
+
 
 def ell_comps_from(axis_ratio: float, angle: float, xp=np) -> Tuple[float, float]:
     """
@@ -72,9 +87,9 @@ def axis_ratio_and_angle_from(
     if xp.__name__.startswith("jax"):
         import jax
 
-        fac = jax.lax.min(fac, 0.999)
+        fac = jax.lax.min(fac, ELL_COMPS_MAGNITUDE_CLAMP)
     else:  # NumPy
-        fac = np.minimum(fac, 0.999)
+        fac = np.minimum(fac, ELL_COMPS_MAGNITUDE_CLAMP)
 
     axis_ratio = (1 - fac) / (1 + fac)
     return axis_ratio, angle

autogalaxy/profiles/geometry_profiles.py

@@ -237,6 +237,33 @@ def __init__(
         validate.validate_ell_comps(ell_comps=ell_comps)
         self.ell_comps = ell_comps
 
+    def __model_constraint__(self, xp=np):
+        """
+        How far beyond the ellipticity clamp this profile's `ell_comps` sit.
+
+        Zero inside the valid region, growing with the distance outside it. This
+        is the traced counterpart of the `validate_ell_comps` guard called in
+        `__init__` above: that guard states the same geometry but signals by
+        raising, which works only for concrete scalars — under a JAX trace the
+        condition is a tracer and a `raise` is impossible, so it returns early
+        (`validate.py:153-154`). PyAutoFit consumes this on the traced path
+        instead, to count multi-start lanes trapped where the clamp has removed
+        the radial gradient.
+
+        Two thresholds are in play and this one is deliberately the clamp's, not
+        the guard's. The guard rejects magnitude >= 1.0, where the geometry stops
+        meaning anything; the clamp saturates at 0.999, where the *gradient*
+        dies. The annulus between them is reachable and is exactly where a
+        gradient search sticks while still passing validation, so keying this to
+        the guard's threshold would miss it.
+
+        Returns a distance rather than a boolean so a future penalty term can use
+        the value directly — it carries a usable gradient back toward the valid
+        region, which a boolean would not.
+        """
+        magnitude = xp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2)
+        return xp.maximum(magnitude - convert.ELL_COMPS_MAGNITUDE_CLAMP, 0.0)
+
     def axis_ratio(self, xp=np) -> float:
         """
         The ratio of the minor-axis to major-axis (b/a) of the ellipse defined by profile (0.0 > q > 1.0).

autogalaxy/profiles/light/standard/sersic.py

@@ -25,6 +25,7 @@
 from autogalaxy.profiles.light.decorators import (
     check_operated_only,
 )
+from autogalaxy import convert
 from autogalaxy.profiles import validate
 
 
@@ -178,7 +179,9 @@ def _eccentric_radii_grid_from_cartesian(
                 1.0e-12,
             )
         )
-        ell_comps_scale = xp.minimum(1.0, 0.999 / ell_comps_norm)
+        ell_comps_scale = xp.minimum(
+            1.0, convert.ELL_COMPS_MAGNITUDE_CLAMP / ell_comps_norm
+        )
 
         ell_comps_y = xp.multiply(ell_comps_y, ell_comps_scale)
         ell_comps_x = xp.multiply(ell_comps_x, ell_comps_scale)

test_autogalaxy/profiles/test_model_constraint.py

@@ -0,0 +1,85 @@
+import numpy as np
+import pytest
+
+import autogalaxy as ag
+from autogalaxy import convert
+from autogalaxy.profiles.geometry_profiles import EllProfile
+
+
+class TestDeclaration:
+    def test_every_elliptical_profile_inherits_it(self):
+        """`ell_comps` has exactly one assignment site, at `EllProfile`, so the
+        declaration reaches every elliptical light and mass profile."""
+        for cls in (ag.mp.Isothermal, ag.mp.PowerLaw, ag.lp.Sersic, ag.lp.Exponential):
+            assert issubclass(cls, EllProfile)
+            assert callable(getattr(cls, "__model_constraint__", None))
+
+    def test_spherical_profiles_inherit_it_and_never_violate(self):
+        """Spherical profiles subclass their elliptical counterpart
+        (`IsothermalSph` -> `Isothermal` -> ... -> `EllProfile`), so they carry
+        the declaration too. Their `ell_comps` are pinned at (0, 0), so it is
+        always satisfied — correct, if a few wasted ops."""
+        profile = ag.mp.IsothermalSph(einstein_radius=1.0)
+        assert profile.ell_comps == (0.0, 0.0)
+        assert profile.__model_constraint__() == 0.0
+
+
+class TestViolation:
+    def test_zero_well_inside(self):
+        profile = ag.mp.Isothermal(ell_comps=(0.3, 0.2))
+        assert profile.__model_constraint__() == 0.0
+
+    def test_zero_just_inside_the_clamp(self):
+        profile = ag.mp.Isothermal(ell_comps=(0.99, 0.0))
+        assert profile.__model_constraint__() == 0.0
+
+    def test_positive_beyond_the_clamp(self):
+        """Constructed inside the guard's valid region (magnitude < 1) but past
+        the clamp — the annulus keyed to 0.999 rather than 1.0 exists for."""
+        profile = ag.mp.Isothermal(ell_comps=(0.9995, 0.0))
+        assert profile.__model_constraint__() > 0.0
+
+    def test_grows_with_distance(self):
+        near = ag.mp.Isothermal(ell_comps=(0.9995, 0.0)).__model_constraint__()
+        far = ag.mp.Isothermal(ell_comps=(0.99999, 0.0)).__model_constraint__()
+        assert far > near
+
+    def test_matches_the_clamp_threshold(self):
+        profile = ag.mp.Isothermal(ell_comps=(0.9995, 0.0))
+        assert profile.__model_constraint__() == pytest.approx(
+            0.9995 - convert.ELL_COMPS_MAGNITUDE_CLAMP
+        )
+
+
+class TestAgreesWithTheGuard:
+    """The constraint and `validate_ell_comps` state the same geometry; the
+    constraint simply engages earlier, at the clamp rather than at validity."""
+
+    def test_guard_rejects_what_the_constraint_flags_beyond_one(self):
+        with pytest.raises(Exception):
+            ag.mp.Isothermal(ell_comps=(1.2, 0.0))
+
+    def test_corner_region_is_flagged_and_rejected(self):
+        """Both components inside (-1, 1), magnitude above 1 — invisible to any
+        per-parameter limit check."""
+        with pytest.raises(Exception):
+            ag.mp.Isothermal(ell_comps=(0.8, 0.8))
+
+    def test_constraint_covers_the_annulus_the_guard_permits(self):
+        """0.999 <= magnitude < 1.0: guard-valid, clamp-saturated."""
+        profile = ag.mp.Isothermal(ell_comps=(0.9995, 0.0))
+        magnitude_squared = 0.9995**2
+        assert magnitude_squared < 1.0  # the guard is satisfied
+        assert profile.__model_constraint__() > 0.0  # the constraint is not
+
+
+class TestClampConstant:
+    def test_conversion_saturates_at_the_constant(self):
+        at = convert.axis_ratio_from(
+            ell_comps=(convert.ELL_COMPS_MAGNITUDE_CLAMP, 0.0)
+        )
+        beyond = convert.axis_ratio_from(ell_comps=(0.99999, 0.0))
+        assert at == pytest.approx(beyond)
+
+    def test_value_is_unchanged_by_the_refactor(self):
+        assert convert.ELL_COMPS_MAGNITUDE_CLAMP == 0.999