diff --git a/autogalaxy/convert.py b/autogalaxy/convert.py index 757b3578..e53d514d 100644 --- a/autogalaxy/convert.py +++ b/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 diff --git a/autogalaxy/profiles/geometry_profiles.py b/autogalaxy/profiles/geometry_profiles.py index c37d32cf..4d2bdd1d 100644 --- a/autogalaxy/profiles/geometry_profiles.py +++ b/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). diff --git a/autogalaxy/profiles/light/standard/sersic.py b/autogalaxy/profiles/light/standard/sersic.py index 8c184b74..4dc4ebc9 100644 --- a/autogalaxy/profiles/light/standard/sersic.py +++ b/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) diff --git a/test_autogalaxy/profiles/test_model_constraint.py b/test_autogalaxy/profiles/test_model_constraint.py new file mode 100644 index 00000000..ce73d53d --- /dev/null +++ b/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