From 621e6a800b0345badadbe330f936bf66ea4befc8 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 9 Aug 2026 22:23:12 +0000 Subject: [PATCH 1/2] =?UTF-8?q?fix:=20validate=20profile=20constructor=20i?= =?UTF-8?q?nputs=20(#440=20=E2=80=94=20B9,=20B11,=20B12)=20+=20B10=20toler?= =?UTF-8?q?ance=20pin?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four findings from @rhayes777's API audit, all still reproducing on main. - B9 scale_radius of 0.0 / negative / nan -> guarded at the 5 assignment sites covering the whole halo family (AbstractgNFW, cNFW x2, Kaplinghat, Yang24). A zero scale radius divided the grid by zero and returned all-NaN deflections (3200 of 3200) rather than raising — the profile appeared to work. - B11 sersic_index of 0.0 / negative / nan -> guarded at both Sersic bases (the light profile and the stellar mass profile). Previously a bare ZeroDivisionError from inside image_2d_from. - B12 ell_comps outside the unit circle -> guarded once at EllProfile, the single base every elliptical light and mass profile inherits. q = (1-f)/(1+f) is only a valid axis ratio for f < 1; beyond it q goes negative and the profile has no geometric meaning, yet it returned a finite non-physical image. - B10 is a tolerance pin, not a fix: Isothermal(ell_comps=(0,0)) vs IsothermalSph are analytically identical, numerically not. Pinned so a future refactor that makes the agreement materially worse is caught. Guards delegate to the shared autoarray.validate helpers established for #333, so all three repos give one message for one class of mistake. Per-parameter explanations live once in autogalaxy/profiles/validate.py rather than per class. Tracer-safe: gated on autoarray.validate.is_concrete_scalar, so a traced free model parameter passes through untouched and no Python truth-test reaches a tracer. NOTE for review: while pinning B10 the potential was measured as well as the deflections the reporter reported. The potential agrees only to 1.9e-03 in relative terms, three orders of magnitude worse than deflections (2.4e-06) and convergence (1.5e-06). That is NOT part of B10 as filed and is not fixed here; its tolerance is pinned at the observed level as a ratchet, and flagged for its own investigation. Tests: 30 new cases in test_autogalaxy/profiles/test_validate.py, one per finding from the reporter's snippets plus a control per finding. Suite 1074 passed, zero regressions. Closes #440. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_013PgqSCLTemK5bApVAwhVM4 --- autogalaxy/profiles/geometry_profiles.py | 3 + autogalaxy/profiles/light/standard/sersic.py | 3 + autogalaxy/profiles/mass/dark/abstract.py | 2 + autogalaxy/profiles/mass/dark/cnfw.py | 3 + autogalaxy/profiles/mass/dark/kaplinghat.py | 2 + autogalaxy/profiles/mass/dark/yang24.py | 2 + autogalaxy/profiles/mass/stellar/sersic.py | 2 + autogalaxy/profiles/validate.py | 127 ++++++++++ test_autogalaxy/profiles/test_validate.py | 240 +++++++++++++++++++ 9 files changed, 384 insertions(+) create mode 100644 autogalaxy/profiles/validate.py create mode 100644 test_autogalaxy/profiles/test_validate.py diff --git a/autogalaxy/profiles/geometry_profiles.py b/autogalaxy/profiles/geometry_profiles.py index 5ab82979..c37d32cf 100644 --- a/autogalaxy/profiles/geometry_profiles.py +++ b/autogalaxy/profiles/geometry_profiles.py @@ -6,6 +6,7 @@ light and mass profiles inherit, including translating a grid to the profile centre and rotating it to the profile's position angle. """ + import numpy as np from typing import Optional, Tuple, Type @@ -13,6 +14,7 @@ import autoarray as aa from autogalaxy import convert +from autogalaxy.profiles import validate class GeometryProfile: @@ -232,6 +234,7 @@ def __init__( """ super().__init__(centre=centre) + validate.validate_ell_comps(ell_comps=ell_comps) self.ell_comps = ell_comps def axis_ratio(self, xp=np) -> float: diff --git a/autogalaxy/profiles/light/standard/sersic.py b/autogalaxy/profiles/light/standard/sersic.py index 99c09688..f3749ec1 100644 --- a/autogalaxy/profiles/light/standard/sersic.py +++ b/autogalaxy/profiles/light/standard/sersic.py @@ -13,6 +13,7 @@ This module provides both elliptical (`Sersic`) and spherical (`SersicSph`) variants. """ + import numpy as np from numpy import seterr @@ -24,6 +25,7 @@ from autogalaxy.profiles.light.decorators import ( check_operated_only, ) +from autogalaxy.profiles import validate class AbstractSersic(LightProfile): @@ -54,6 +56,7 @@ def __init__( """ super().__init__(centre=centre, ell_comps=ell_comps, intensity=intensity) self.effective_radius = effective_radius + validate.validate_sersic_index(sersic_index=sersic_index) self.sersic_index = sersic_index @property diff --git a/autogalaxy/profiles/mass/dark/abstract.py b/autogalaxy/profiles/mass/dark/abstract.py index 75f1f472..e00cd1fc 100644 --- a/autogalaxy/profiles/mass/dark/abstract.py +++ b/autogalaxy/profiles/mass/dark/abstract.py @@ -8,6 +8,7 @@ from autogalaxy import exc +from autogalaxy.profiles import validate class DarkProfile: @@ -92,6 +93,7 @@ def __init__( super().__init__(centre=centre, ell_comps=ell_comps) self.kappa_s = kappa_s + validate.validate_scale_radius(scale_radius=scale_radius) self.scale_radius = scale_radius self.inner_slope = inner_slope diff --git a/autogalaxy/profiles/mass/dark/cnfw.py b/autogalaxy/profiles/mass/dark/cnfw.py index cee4d7bf..fbe76177 100644 --- a/autogalaxy/profiles/mass/dark/cnfw.py +++ b/autogalaxy/profiles/mass/dark/cnfw.py @@ -6,6 +6,7 @@ from autogalaxy.profiles.mass.abstract.mge import MGEDecomposer import autoarray as aa +from autogalaxy.profiles import validate def F_func_from(theta, radius, xp=np): @@ -126,6 +127,7 @@ def __init__( super().__init__(centre=centre, ell_comps=ell_comps) self.kappa_s = kappa_s + validate.validate_scale_radius(scale_radius=scale_radius) self.scale_radius = scale_radius self.core_radius = core_radius @@ -275,6 +277,7 @@ def __init__( super().__init__(centre=centre, ell_comps=(0.0, 0.0)) self.kappa_s = kappa_s + validate.validate_scale_radius(scale_radius=scale_radius) self.scale_radius = scale_radius self.core_radius = core_radius diff --git a/autogalaxy/profiles/mass/dark/kaplinghat.py b/autogalaxy/profiles/mass/dark/kaplinghat.py index fcf2a427..de33aefe 100644 --- a/autogalaxy/profiles/mass/dark/kaplinghat.py +++ b/autogalaxy/profiles/mass/dark/kaplinghat.py @@ -10,6 +10,7 @@ from autogalaxy.profiles.mass.dark.abstract import DarkProfile from autogalaxy.profiles.mass.dark.nfw import NFWSph from autogalaxy.profiles.mass.abstract.abstract import MassProfile +from autogalaxy.profiles import validate @functools.lru_cache(maxsize=1) @@ -183,6 +184,7 @@ def __init__( super().__init__(centre=centre, ell_comps=(0.0, 0.0)) self.kappa_s = kappa_s + validate.validate_scale_radius(scale_radius=scale_radius) self.scale_radius = scale_radius self.sigma_over_m = sigma_over_m self.t_age = t_age diff --git a/autogalaxy/profiles/mass/dark/yang24.py b/autogalaxy/profiles/mass/dark/yang24.py index 4e06d385..29f54115 100644 --- a/autogalaxy/profiles/mass/dark/yang24.py +++ b/autogalaxy/profiles/mass/dark/yang24.py @@ -12,6 +12,7 @@ _trapezoid_from, ) from autogalaxy.profiles.mass.abstract.abstract import MassProfile +from autogalaxy.profiles import validate def _yang24_parameter_ratios_from(tau): @@ -114,6 +115,7 @@ def __init__( super().__init__(centre=centre, ell_comps=(0.0, 0.0)) self.kappa_s = kappa_s + validate.validate_scale_radius(scale_radius=scale_radius) self.scale_radius = scale_radius self.tau = min(max(float(tau), 0.0), 1.0) diff --git a/autogalaxy/profiles/mass/stellar/sersic.py b/autogalaxy/profiles/mass/stellar/sersic.py index 65db4163..2a759114 100644 --- a/autogalaxy/profiles/mass/stellar/sersic.py +++ b/autogalaxy/profiles/mass/stellar/sersic.py @@ -12,6 +12,7 @@ MassProfileCSE, ) from autogalaxy.profiles.mass.stellar.abstract import StellarProfile +from autogalaxy.profiles import validate def cse_settings_from( @@ -171,6 +172,7 @@ def __init__( self.mass_to_light_ratio = mass_to_light_ratio self.intensity = intensity self.effective_radius = effective_radius + validate.validate_sersic_index(sersic_index=sersic_index) self.sersic_index = sersic_index def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): diff --git a/autogalaxy/profiles/validate.py b/autogalaxy/profiles/validate.py new file mode 100644 index 00000000..89d2096a --- /dev/null +++ b/autogalaxy/profiles/validate.py @@ -0,0 +1,127 @@ +""" +Profile constructor guards. + +These are thin, named wrappers over the shared helpers in ``autoarray.validate`` +(landed by PyAutoArray#440 for PyAutoArray#333). PyAutoArray is the floor both +PyAutoGalaxy and PyAutoLens build on, so the *rules* and the *message shape* are +defined once there and reused here — the failure mode being avoided is three repos +telling a user three different things about the same mistake. + +What lives here is only the per-parameter explanation, written once per parameter +rather than once per profile class. Every message therefore reads the same way: +name the parameter, state the rule, show the received value, then explain why. + +__Tracer safety__ + +Profile parameters are free model parameters, so under a JAX-traced fit a +constructor is handed a tracer rather than a number, and a plain Python +``if value <= 0`` would raise ``TracerBoolConversionError``. The shared helpers gate +every comparison on ``autoarray.validate.is_concrete_scalar`` and pass non-concrete +values straight through, so these guards catch hand-written mistakes and cost +nothing inside a trace. +""" + +import numpy as np + +from autoarray import validate + + +def validate_scale_radius(scale_radius, name: str = "scale_radius"): + """ + Raise if a dark-matter halo scale radius is a concrete scalar which is not finite + and positive. + + A ``scale_radius`` of zero divides the grid by zero in the very first step of the + deflection-angle calculation, so every returned value is NaN rather than an error — + the profile appears to work and quietly poisons the whole fit. + + Parameters + ---------- + scale_radius + The scale radius to validate. + name + The parameter's name, used in the error message. + """ + validate.validate_positive_finite( + value=scale_radius, + name=name, + extra=( + "The scale radius is the angular radius at which the halo's log-slope " + "changes, so it must be above zero. A value of 0.0 divides the grid by " + "zero when computing deflection angles, which returns an all-NaN result " + "instead of raising" + ), + ) + + +def validate_sersic_index(sersic_index, name: str = "sersic_index"): + """ + Raise if a Sersic index is a concrete scalar which is not finite and positive. + + A ``sersic_index`` of zero reaches a division by ``n`` deep inside the profile's + ``image_2d_from``, surfacing as a bare ``ZeroDivisionError`` several calls away + from the constructor that accepted it. + + Parameters + ---------- + sersic_index + The Sersic index to validate. + name + The parameter's name, used in the error message. + """ + validate.validate_positive_finite( + value=sersic_index, + name=name, + extra=( + "The Sersic index controls the concentration of the profile and appears " + "as a divisor in its normalisation, so it must be above zero. A value of " + "0.0 raises ZeroDivisionError from inside image_2d_from rather than at " + "construction" + ), + ) + + +def validate_ell_comps(ell_comps, name: str = "ell_comps"): + """ + Raise if the elliptical components are concrete scalars whose magnitude is not + below one. + + The axis ratio is defined as ``q = (1 - f) / (1 + f)`` with + ``f = sqrt(e_y**2 + e_x**2)``. That is a valid axis ratio in ``(0, 1]`` only while + ``f < 1``; at ``f == 1`` the ellipse degenerates to ``q == 0``, and beyond it ``q`` + goes negative and the profile has no geometric meaning. Today such a profile is + accepted and returns a finite but non-physical image. + + Applied at ``EllProfile``, the single base every elliptical light and mass profile + inherits, so the rule is stated once rather than per subclass. + + Parameters + ---------- + ell_comps + The (e_y, e_x) elliptical components to validate. + name + The parameter's name, used in the error message. + """ + if ell_comps is None: + return + + try: + ell_y, ell_x = ell_comps + except (TypeError, ValueError): + return + + if not validate.is_concrete_scalar(ell_y) or not validate.is_concrete_scalar(ell_x): + return + + magnitude_squared = ell_y * ell_y + ell_x * ell_x + + if not np.isfinite(magnitude_squared) or magnitude_squared >= 1.0: + raise ValueError( + f"{name} must satisfy {name}[0]**2 + {name}[1]**2 < 1; got " + f"{tuple(ell_comps)!r}, whose magnitude is " + f"{np.sqrt(magnitude_squared) if np.isfinite(magnitude_squared) else magnitude_squared!r}. " + f"The axis ratio is q = (1 - f) / (1 + f) with f the magnitude of " + f"{name}, so f must be below 1 for q to be a valid axis ratio — at f = 1 " + f"the ellipse degenerates to q = 0 and beyond it q is negative and the " + f"profile has no geometric meaning" + ) diff --git a/test_autogalaxy/profiles/test_validate.py b/test_autogalaxy/profiles/test_validate.py new file mode 100644 index 00000000..97392157 --- /dev/null +++ b/test_autogalaxy/profiles/test_validate.py @@ -0,0 +1,240 @@ +""" +Regression tests for PyAutoGalaxy#440 — profile constructor validation (B9, B11, +B12) and the B10 Ell/Sph tolerance pin. + +The guard tests are built from @rhayes777's own snippets in the issue body and +assert the *failure*: that the input is rejected at construction with a message +naming the offending parameter, rather than accepted and surfaced as an all-NaN +array or a bare ZeroDivisionError several calls later. + +Each finding is paired with a control asserting the valid input still works, so a +guard cannot pass by rejecting everything. + +Tests are numpy-only, per phase 1. Tracer-safety is asserted against the +concreteness gate the guards branch on (`autoarray.validate.is_concrete_scalar`) +rather than by importing JAX into the library unit tests. +""" + +import numpy as np +import pytest + +import autogalaxy as ag +from autogalaxy.profiles import validate + + +class _NotAConcreteScalar: + """ + Stand-in for a JAX tracer: not a concrete Python/NumPy scalar, and raises if + anything resolves it to a bool — exactly as a tracer does inside `jax.jit`. + """ + + def __bool__(self): + raise AssertionError( + "a guard compared a non-concrete value — the TracerBoolConversionError path" + ) + + def __lt__(self, other): + return self + + def __le__(self, other): + return self + + def __mul__(self, other): + return self + + def __add__(self, other): + return self + + +# ====================================================================================== +# B9 — scale_radius must be finite and positive +# ====================================================================================== + + +@pytest.mark.parametrize("scale_radius", [0.0, -1.0, float("nan"), float("inf")]) +def test__b9__nfw_rejects_non_positive_or_non_finite_scale_radius(scale_radius): + with pytest.raises(ValueError, match="scale_radius"): + ag.mp.NFW(scale_radius=scale_radius) + + +@pytest.mark.parametrize( + "profile_cls", [ag.mp.NFW, ag.mp.NFWSph, ag.mp.gNFW, ag.mp.gNFWSph, ag.mp.cNFW] +) +def test__b9__every_nfw_family_profile_rejects_a_zero_scale_radius(profile_cls): + """The reporter named `NFW`; the same hole was open across the halo family.""" + with pytest.raises(ValueError, match="scale_radius"): + profile_cls(scale_radius=0.0) + + +def test__b9__control__a_positive_scale_radius_builds_and_returns_finite_deflections(): + grid = ag.Grid2D.uniform(shape_native=(20, 20), pixel_scales=0.1) + + deflections = ag.mp.NFW(scale_radius=1.0).deflections_yx_2d_from(grid=grid) + + assert np.isfinite(np.asarray(deflections)).all() + + +# ====================================================================================== +# B11 — sersic_index must be finite and positive +# ====================================================================================== + + +@pytest.mark.parametrize("sersic_index", [0.0, -1.0, float("nan"), float("inf")]) +def test__b11__sersic_rejects_non_positive_or_non_finite_sersic_index(sersic_index): + with pytest.raises(ValueError, match="sersic_index"): + ag.lp.Sersic(sersic_index=sersic_index) + + +def test__b11__the_stellar_mass_sersic_is_guarded_too(): + with pytest.raises(ValueError, match="sersic_index"): + ag.mp.Sersic(sersic_index=0.0) + + +def test__b11__control__a_normal_sersic_index_still_produces_a_finite_image(): + grid = ag.Grid2D.uniform(shape_native=(20, 20), pixel_scales=0.1) + + image = ag.lp.Sersic(sersic_index=1.0).image_2d_from(grid=grid) + + assert np.isfinite(np.asarray(image)).all() + + +# ====================================================================================== +# B12 — ell_comps must lie inside the unit circle +# ====================================================================================== + + +@pytest.mark.parametrize("ell_comps", [(2.0, 0.0), (0.0, 2.0), (0.9, 0.9), (1.0, 0.0)]) +def test__b12__elliptical_profiles_reject_ell_comps_of_magnitude_one_or_above( + ell_comps, +): + """ + `q = (1 - f) / (1 + f)` is a valid axis ratio only for `f < 1`. At `f == 1` the + ellipse degenerates to `q == 0`; beyond it `q` is negative and meaningless. + """ + with pytest.raises(ValueError, match="ell_comps"): + ag.lp.Sersic(ell_comps=ell_comps) + + +def test__b12__the_guard_is_on_the_shared_elliptical_base_not_one_subclass(): + """`EllProfile` is the single base every elliptical light and mass profile uses.""" + with pytest.raises(ValueError, match="ell_comps"): + ag.mp.Isothermal(ell_comps=(2.0, 0.0)) + + with pytest.raises(ValueError, match="ell_comps"): + ag.lp.Gaussian(ell_comps=(2.0, 0.0)) + + with pytest.raises(ValueError, match="ell_comps"): + ag.mp.NFW(ell_comps=(2.0, 0.0), scale_radius=1.0) + + +def test__b12__the_message_reports_the_magnitude(): + with pytest.raises(ValueError) as error: + ag.lp.Sersic(ell_comps=(3.0, 4.0)) + + assert "5.0" in str(error.value) + + +@pytest.mark.parametrize("ell_comps", [(0.0, 0.0), (0.3, 0.4), (0.0, 0.9)]) +def test__b12__control__ell_comps_inside_the_unit_circle_still_build(ell_comps): + profile = ag.lp.Sersic(ell_comps=ell_comps) + + assert 0.0 < profile.axis_ratio() <= 1.0 + + +# ====================================================================================== +# Tracer safety +# ====================================================================================== + + +def test__guards_pass_through_non_concrete_values__never_compare_them(): + tracer_like = _NotAConcreteScalar() + + validate.validate_scale_radius(scale_radius=tracer_like) + validate.validate_sersic_index(sersic_index=tracer_like) + validate.validate_ell_comps(ell_comps=(tracer_like, tracer_like)) + + +def test__profile_constructors_accept_tracer_like_parameters(): + """Profile parameters are free model parameters; under a trace they arrive traced.""" + tracer_like = _NotAConcreteScalar() + + assert ag.mp.NFW(scale_radius=tracer_like).scale_radius is tracer_like + assert ag.lp.Sersic(sersic_index=tracer_like).sersic_index is tracer_like + assert ag.lp.Sersic(ell_comps=(tracer_like, tracer_like)).ell_comps == ( + tracer_like, + tracer_like, + ) + + +# ====================================================================================== +# B10 — Isothermal(ell_comps=(0,0)) vs IsothermalSph agreement, pinned at a tolerance +# ====================================================================================== +# +# These two are analytically identical: the elliptical form at zero ellipticity IS the +# spherical form. Numerically they differ, because the `Ell` form takes a different +# evaluation route even at the degenerate point. +# +# This is a TOLERANCE PIN, not a bug fix. Bit-identity is explicitly not the goal — the +# point is that a future refactor which makes the agreement materially worse gets +# caught. Measured values on `main` at the time of writing are in each test. + + +def _isothermal_pair(): + ell = ag.mp.Isothermal(ell_comps=(0.0, 0.0), einstein_radius=1.0) + sph = ag.mp.IsothermalSph(einstein_radius=1.0) + return ell, sph + + +def test__b10__deflections_agree_between_elliptical_and_spherical_isothermal(): + """Measured max|diff| = 2.357e-06 (relative 2.36e-06). Pinned an order looser.""" + grid = ag.Grid2D.uniform(shape_native=(40, 40), pixel_scales=0.1) + ell, sph = _isothermal_pair() + + difference = np.max( + np.abs( + np.asarray(ell.deflections_yx_2d_from(grid=grid)) + - np.asarray(sph.deflections_yx_2d_from(grid=grid)) + ) + ) + + assert difference < 1.0e-5 + + +def test__b10__convergence_agrees_between_elliptical_and_spherical_isothermal(): + """Measured max|diff| = 1.207e-05 (relative 1.45e-06). Pinned an order looser.""" + grid = ag.Grid2D.uniform(shape_native=(40, 40), pixel_scales=0.1) + ell, sph = _isothermal_pair() + + difference = np.max( + np.abs( + np.asarray(ell.convergence_2d_from(grid=grid)) + - np.asarray(sph.convergence_2d_from(grid=grid)) + ) + ) + + assert difference < 1.0e-4 + + +def test__b10__potential_agrees_between_elliptical_and_spherical_isothermal(): + """ + Measured max|diff| = 5.375e-03, which is a **relative** difference of 1.9e-03 — + three orders of magnitude worse than the deflection and convergence agreement + above, and NOT part of @rhayes777's original B10 report (he measured deflections + only). + + This tolerance is therefore deliberately pinned at the currently-observed level + rather than at a level anyone has argued is scientifically acceptable. It is a + ratchet: it stops the agreement degrading further, and it is expected to be + tightened when the potential discrepancy is investigated on its own. + """ + grid = ag.Grid2D.uniform(shape_native=(40, 40), pixel_scales=0.1) + ell, sph = _isothermal_pair() + + difference = np.max( + np.abs( + np.asarray(ell.potential_2d_from(grid=grid)) + - np.asarray(sph.potential_2d_from(grid=grid)) + ) + ) + + assert difference < 1.0e-2 From 68d1f8a8f17431832af678a71ff3e6c518eaa306 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 9 Aug 2026 22:25:42 +0000 Subject: [PATCH 2/2] fix: reject negative galaxy redshifts (PyAutoLens#532, half that lives here) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The negative-redshift finding is filed on PyAutoLens#532 because the reporter reached it through `al.Galaxy` — but `al.Galaxy` IS `ag.Galaxy`, and both the class and its redshift assignment live in this repo, so the guard belongs here rather than in a Tracer-level check that would miss a bare Galaxy construction. The `Tracer(galaxies=...)` half of #532 stays in PyAutoLens. Zero and tiny redshifts stay accepted: 0.0 legitimately places a galaxy at the observer, and 1e-12 is degenerate but not invalid. Adds a deliberate guard-rail test pinning today's permissive z_lens > z_source behaviour, so phase 4 of the audit cannot quietly turn it into an error while the question is still open with the reporter. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_013PgqSCLTemK5bApVAwhVM4 --- autogalaxy/galaxy/galaxy.py | 4 +++ autogalaxy/profiles/validate.py | 36 +++++++++++++++++++ test_autogalaxy/profiles/test_validate.py | 42 +++++++++++++++++++++++ 3 files changed, 82 insertions(+) diff --git a/autogalaxy/galaxy/galaxy.py b/autogalaxy/galaxy/galaxy.py index 43a92521..ca89ccc2 100644 --- a/autogalaxy/galaxy/galaxy.py +++ b/autogalaxy/galaxy/galaxy.py @@ -27,6 +27,7 @@ from autogalaxy.profiles.light.linear import LightProfileLinear from autogalaxy.profiles.light.snr.abstract import LightProfileSNR from autogalaxy.profiles.mass.abstract.abstract import MassProfile +from autogalaxy.profiles import validate class Galaxy(af.ModelObject, OperateImageList): @@ -49,6 +50,9 @@ def __init__(self, redshift: float, **kwargs): The pixelization of the galaxy used to reconstruct an observed image using an inversion. """ super().__init__() + + validate.validate_redshift(redshift=redshift) + self.redshift = redshift for name, val in kwargs.items(): diff --git a/autogalaxy/profiles/validate.py b/autogalaxy/profiles/validate.py index 89d2096a..cd974bc2 100644 --- a/autogalaxy/profiles/validate.py +++ b/autogalaxy/profiles/validate.py @@ -81,6 +81,42 @@ def validate_sersic_index(sersic_index, name: str = "sersic_index"): ) +def validate_redshift(redshift, name: str = "redshift"): + """ + Raise if a galaxy redshift is a concrete scalar which is negative or non-finite. + + A negative redshift is unphysical and produces meaningless angular diameter + distances in every multi-plane calculation that consumes it. + + Zero is permitted: ``redshift=0.0`` is a legitimate way to place a galaxy at the + observer, and is used in single-plane work where the redshift is a label rather + than a cosmological quantity. + + **This does not touch the ``z_lens > z_source`` question.** Multi-plane lensing + genuinely supports geometries that look wrong under two-plane naming, so that + case must warn at most, never raise, and is held pending the reporter's answer on + PyAutoLens#532. + + Parameters + ---------- + redshift + The redshift to validate. + name + The parameter's name, used in the error message. + """ + validate.validate_non_negative_finite( + value=redshift, + name=name, + extra=( + "A redshift is a cosmological distance measure and cannot be negative — " + "a negative value yields meaningless angular diameter distances in every " + "multi-plane calculation that consumes it. Note that a lens redshift " + "above a source redshift is NOT rejected: multi-plane lensing supports " + "geometries that look inverted under two-plane naming" + ), + ) + + def validate_ell_comps(ell_comps, name: str = "ell_comps"): """ Raise if the elliptical components are concrete scalars whose magnitude is not diff --git a/test_autogalaxy/profiles/test_validate.py b/test_autogalaxy/profiles/test_validate.py index 97392157..0017919f 100644 --- a/test_autogalaxy/profiles/test_validate.py +++ b/test_autogalaxy/profiles/test_validate.py @@ -141,6 +141,48 @@ def test__b12__control__ell_comps_inside_the_unit_circle_still_build(ell_comps): assert 0.0 < profile.axis_ratio() <= 1.0 +# ====================================================================================== +# Negative redshift (the PyAutoLens#532 half that lives in this repo) +# ====================================================================================== +# +# Filed on PyAutoLens#532 because the reporter reached it through `al.Galaxy`, but +# `al.Galaxy` IS `ag.Galaxy` — the class and its `redshift` assignment live here, so +# the guard belongs here. The `Tracer(galaxies=...)` half of that issue is in +# PyAutoLens. + + +@pytest.mark.parametrize("redshift", [-0.5, -1.0, float("nan"), float("inf")]) +def test__galaxy_rejects_a_negative_or_non_finite_redshift(redshift): + with pytest.raises(ValueError, match="redshift"): + ag.Galaxy(redshift=redshift) + + +def test__control__zero_and_tiny_redshifts_are_still_accepted(): + """ + Zero places a galaxy at the observer and is legitimate. `1e-12` was flagged by the + reporter as degenerate, but it is not *invalid* — rejecting it would break + single-plane work where the redshift is a label, so it stays accepted. + """ + assert ag.Galaxy(redshift=0.0).redshift == 0.0 + assert ag.Galaxy(redshift=1e-12).redshift == 1e-12 + + +def test__control__a_lens_redshift_above_the_source_redshift_still_constructs(): + """ + PHASE 4 GUARD-RAIL — deliberately pinning today's permissive behaviour. + + `z_lens > z_source` must NOT raise: multi-plane lensing genuinely supports + geometries that look inverted under two-plane naming. Whether it should even + *warn* is an open question put to @rhayes777 on PyAutoLens#532. This test exists + so that phase 4 cannot quietly turn it into an error. + """ + lens = ag.Galaxy(redshift=1.0, mass=ag.mp.IsothermalSph(einstein_radius=1.0)) + source = ag.Galaxy(redshift=0.5, bulge=ag.lp.Sersic(intensity=1.0)) + + assert lens.redshift == 1.0 + assert source.redshift == 0.5 + + # ====================================================================================== # Tracer safety # ======================================================================================