From de3b43663e0977cb3ba5a04f687d9949c7115e58 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 28 Jul 2026 16:57:20 +0100 Subject: [PATCH 1/2] fix: PointSolver crashes on zero images and on too-coarse precision Two realistic inputs crashed `PointSolver.solve` with errors naming nothing the caller controls. 1. A source-plane coordinate outside the region the image-plane grid tiles produces no images. `filtered_means` is then empty and the `Grid2DIrregular` built from it exposes `.array` as a bare list, so the `axis=1` reductions raised `numpy.exceptions.AxisError: axis 1 is out of bounds for array of dimension 1`. Zero images is a legitimate answer, so the solver now returns a correctly shaped empty grid and logs why. (The report expected a single image here; the tiling genuinely finds none -- the coordinate lies outside the tiled region.) 2. `n_steps` is `ceil(log2(scale / pixel_scale_precision))`, which goes negative once the requested precision is coarser than the initial triangle scale. The existing guard tested `== 0`, so negative values slipped through, `range(-1)` yielded no steps, and `steps[-1]` raised `IndexError: list index out of range`. The guard now tests `<= 0` and reports `pixel_scale_precision`, the triangle scale, and a workable value. It deliberately does not clamp to 1 step, which would silently solve at a precision the caller did not ask for. Reported by @rhayes777 in #531. Phase 1 of PyAutoArray#416 (epic PyAutoArray#415). Co-Authored-By: Claude Opus 5 --- autolens/point/solver/point_solver.py | 21 +++ autolens/point/solver/shape_solver.py | 22 +++- .../point/triangles/test_solver_edge_cases.py | 121 ++++++++++++++++++ 3 files changed, 162 insertions(+), 2 deletions(-) create mode 100644 test_autolens/point/triangles/test_solver_edge_cases.py diff --git a/autolens/point/solver/point_solver.py b/autolens/point/solver/point_solver.py index dfb3e393c..8c28cf79e 100644 --- a/autolens/point/solver/point_solver.py +++ b/autolens/point/solver/point_solver.py @@ -122,6 +122,27 @@ def solve( tracer=tracer, points=kept_triangles.means, xp=xp ) + # When no triangle traces to the source-plane coordinate -- e.g. the coordinate lies + # outside the region the image-plane grid tiles -- `filtered_means` is empty. The + # `Grid2DIrregular` built from it exposes `.array` as a bare list rather than a 2D array, + # so the `axis=1` reductions below used to raise + # `numpy.exceptions.AxisError: axis 1 is out of bounds for array of dimension 1`. + # + # Zero images is a legitimate answer, not a failure, so return a correctly-shaped empty + # grid. `filtered_means` is a Python list, so its length is static (including under a + # `jax.jit` trace) and this branch is trace-safe. The JAX path pads to + # `MAX_CONTAINING_SIZE` and never reaches it. + if len(filtered_means) == 0: + + logger.warning( + f"PointSolver.solve found no images for source-plane coordinate " + f"{tuple(source_plane_coordinate)}, so an empty grid is returned. This usually " + f"means the coordinate lies outside the region traced by the image-plane grid, or " + f"that every candidate image was removed by the magnification threshold." + ) + + return aa.Grid2DIrregular(xp.zeros((0, 2))) + solution = aa.Grid2DIrregular( [pair for pair in filtered_means], xp=xp ).array diff --git a/autolens/point/solver/shape_solver.py b/autolens/point/solver/shape_solver.py index 7b155b3dc..161db3bf4 100644 --- a/autolens/point/solver/shape_solver.py +++ b/autolens/point/solver/shape_solver.py @@ -285,9 +285,27 @@ def solve_triangles( ------- A list of image plane coordinates that are traced to the source plane coordinate. """ - if self.n_steps == 0: + # `n_steps` is `ceil(log2(scale / pixel_scale_precision))`, which is <= 0 whenever the + # requested precision is coarser than (or equal to) the initial triangle scale. A negative + # value used to slip past an `== 0` check and make `steps` an empty list, surfacing as + # `IndexError: list index out of range` on `steps[-1]` below. Reject it here, naming the + # parameter the caller controls -- do NOT clamp to 1, which would silently solve at a + # precision the caller did not ask for. + if self.n_steps <= 0: raise ValueError( - "The target pixel scale is too large to subdivide the triangles." + f""" + The requested `pixel_scale_precision` is too large to subdivide the triangles. + + pixel_scale_precision = {self.pixel_scale_precision} + initial triangle scale = {self.scale} + + The solver refines triangles by repeated bisection, so it needs + `pixel_scale_precision` to be smaller than the initial triangle scale; here it + would require {self.n_steps} subdivision steps. + + Decrease `pixel_scale_precision` (e.g. to {self.scale / 10.0:.3g} or smaller) so the + solver can resolve an image. + """ ) steps = list( diff --git a/test_autolens/point/triangles/test_solver_edge_cases.py b/test_autolens/point/triangles/test_solver_edge_cases.py new file mode 100644 index 000000000..0dc2f8fd7 --- /dev/null +++ b/test_autolens/point/triangles/test_solver_edge_cases.py @@ -0,0 +1,121 @@ +""" +Regression tests for @rhayes777's audit finding in PyAutoLens#531. + +Two realistic inputs crashed ``PointSolver.solve`` with errors that named nothing the caller +controls: + +- a source-plane coordinate outside the region the image-plane grid tiles + -> ``numpy.exceptions.AxisError: axis 1 is out of bounds for array of dimension 1`` +- a ``pixel_scale_precision`` coarser than the initial triangle scale + -> ``IndexError: list index out of range`` + +Both are inputs a user reaches on purpose: configurations outside the caustic are routine during +model exploration, and loosening the precision for a quick first pass is an obvious thing to +script. +""" + +import numpy as np +import pytest + +import autolens as al + + +@pytest.fixture +def solver_grid(): + return al.Grid2D.uniform(shape_native=(80, 80), pixel_scales=0.05) + + +@pytest.fixture +def lens_galaxy(): + return al.Galaxy( + redshift=0.5, + mass=al.mp.Isothermal(centre=(0.0, 0.0), ell_comps=(0.1, 0.0), einstein_radius=1.0), + ) + + +def _tracer(lens_galaxy, source_plane_coordinate): + source = al.Galaxy( + redshift=1.0, point_0=al.ps.Point(centre=source_plane_coordinate) + ) + return al.Tracer(galaxies=[lens_galaxy, source]) + + +def test__source_outside_tiled_region__returns_empty_grid(solver_grid, lens_galaxy): + """ + No image is a legitimate answer, so the solver returns a correctly-shaped empty grid rather + than raising `AxisError` from the `axis=1` reductions. + """ + + source_plane_coordinate = (5.0, 5.0) + + solver = al.PointSolver.for_grid(grid=solver_grid, pixel_scale_precision=0.001) + + result = solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + assert len(result) == 0 + assert np.asarray(result.array).shape == (0, 2) + + +def test__source_outside_tiled_region__warns(solver_grid, lens_galaxy, caplog): + """The empty result is explained, so it does not read as a silent no-op.""" + + source_plane_coordinate = (5.0, 5.0) + + solver = al.PointSolver.for_grid(grid=solver_grid, pixel_scale_precision=0.001) + + with caplog.at_level("WARNING"): + solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + assert "found no images" in caplog.text + + +@pytest.mark.parametrize("pixel_scale_precision", [0.05, 0.1, 0.2, 0.5]) +def test__precision_coarser_than_triangle_scale__raises_naming_the_parameter( + solver_grid, lens_galaxy, pixel_scale_precision +): + """ + `n_steps` is `ceil(log2(scale / pixel_scale_precision))`, which is <= 0 once the precision + reaches the triangle scale. A negative value used to slip past an `== 0` check and produce + `IndexError: list index out of range`. + + The grid's `pixel_scales` is 0.05, so every value here is at or beyond the boundary. + """ + + source_plane_coordinate = (0.05, 0.02) + + solver = al.PointSolver.for_grid( + grid=solver_grid, pixel_scale_precision=pixel_scale_precision + ) + + with pytest.raises(ValueError) as error: + solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + # the message must name the parameter the caller controls + assert "pixel_scale_precision" in str(error.value) + + +def test__precision_fine_enough__still_solves(solver_grid, lens_galaxy): + """ + The control from the report: the same source at a workable precision returns its four images. + Guards the fix against over-reach. + """ + + source_plane_coordinate = (0.05, 0.02) + + solver = al.PointSolver.for_grid(grid=solver_grid, pixel_scale_precision=0.001) + + result = solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + assert len(result) == 4 From 37bb880c80173bb88723ac662895619b8cd14884 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Tue, 28 Jul 2026 17:24:23 +0100 Subject: [PATCH 2/2] fix: give FitPositionsImagePairAll a no-image floor; recover from zero images Follow-up to de3b4366, from an independent Codex review of that commit. Making `PointSolver.solve` return an empty grid rather than raising made two pre-existing consumer gaps reachable; both are fixed here. 1. `FitPositionsImagePairAll.chi_squared` returned NaN with no model positions: `n_permutations` is 0, so `-log(0)` is +inf while the permutation sum is -inf. `fitness.py` converts a NaN log-likelihood into `resample_figure_of_merit`, so the model was silently resampled instead of scored -- precisely what `no_image_residual` exists to prevent. `FitPositionsImagePair` and `FitPositionsImagePairRepeat` both already applied that floor; `PairAll` was the one sibling that did not. It now returns the same value they do (verified equal, not merely finite). Selected with `xp.where`, not a Python `if`: the model-position count is a traced value under `jax.jit`. The log is taken on a clamped count so the NaN never forms in the branch `where` discards. 2. `Result.image_plane_multiple_image_positions` invoked its inward-walk recovery only for exactly one image, so zero images fell through as an empty grid and `SourceMaxSeparation` reduced over it with `max()` on an empty sequence. Zero is the case that most needs the recovery; `== 1` -> `<= 1`. Also from the same review: the empty-result warning now fires on the final result, catching the second route to empty (all candidates rejected by `magnification_threshold`, where `_filter_low_magnification` preserves the array length and writes NaN rows). Corrects a comment that described `filtered_means` as a Python list -- it is an array -- and makes the `AbstractSolver` guard message shape-agnostic, since `ShapeSolver` shares it. Co-Authored-By: Claude Opus 5 --- autolens/analysis/result.py | 9 ++- .../point/fit/positions/image/pair_all.py | 35 ++++++++++- autolens/point/solver/point_solver.py | 62 ++++++++++++------- autolens/point/solver/shape_solver.py | 2 +- .../fit/positions/image/test_pair_all.py | 57 +++++++++++++++++ .../point/triangles/test_solver_edge_cases.py | 47 ++++++++++++++ 6 files changed, 187 insertions(+), 25 deletions(-) diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index 8f4efb5b6..6f3c81810 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -119,7 +119,14 @@ def image_plane_multiple_image_positions( plane_redshift=plane_redshift, ) - if multiple_images.shape[0] == 1: + # `<= 1` rather than `== 1`: the solver can legitimately return zero images (the + # source-plane coordinate falls outside the region the image-plane grid tiles, or every + # candidate is rejected by the magnification threshold), and that is the case that most + # needs the inward-walk recovery below. Zero used to be unreachable here because + # `PointSolver.solve` raised before returning; now it returns an empty grid, and letting + # it through would give `SourceMaxSeparation` an empty array to reduce over + # (`max()` on an empty sequence). + if multiple_images.shape[0] <= 1: return self.image_plane_multiple_image_positions_for_single_image_from() return aa.Grid2DIrregular(values=multiple_images) diff --git a/autolens/point/fit/positions/image/pair_all.py b/autolens/point/fit/positions/image/pair_all.py index 59c255ca2..dbb651db3 100644 --- a/autolens/point/fit/positions/image/pair_all.py +++ b/autolens/point/fit/positions/image/pair_all.py @@ -61,6 +61,11 @@ class FitPositionsImagePairAll(AbstractFitPositionsImagePair): tracer via name pairing if that profile is not found. """ + # The floor every observed position contributes when the solver returns no images, matching + # `FitPositionsImagePair` and `FitPositionsImagePairRepeat`: loudly bad, but finite, so the + # model is scored rather than silently resampled. + no_image_residual = 1.0e4 + def log_p( self, data_position: np.ndarray, @@ -150,12 +155,38 @@ def chi_squared(self) -> float: self.model_data.array, ).any(axis=1) ) - n_permutations = n_non_nan_model_positions ** len(self.data) - return -2.0 * ( + + # With no finite model positions `n_permutations` is 0, so `-log(0)` is `+inf` while the + # permutation sum is `-inf`, and the two combine to NaN. `fitness.py` converts a NaN + # log-likelihood into `resample_figure_of_merit`, so the model would be silently rejected + # rather than scored -- the exact outcome `no_image_residual` exists to avoid. + # + # `FitPositionsImagePair` and `FitPositionsImagePairRepeat` both fall back to the + # `no_image_residual` floor for every observed position; do the same here, on the same + # (residual / noise) ** 2 scale their chi-squared ends up on. + # + # `n_non_nan_model_positions` is a traced value under `jax.jit`, so this selects with + # `xp.where` rather than a Python `if`. `where` evaluates both branches, so the log is + # taken on a clamped count to keep the NaN from forming in the discarded branch. + noise_map = self._xp.asarray(np.asarray(self.noise_map)) + + no_image_chi_squared = self._xp.sum( + (self.no_image_residual / noise_map) ** 2.0 + ) + + has_image = n_non_nan_model_positions > 0 + + n_permutations = ( + self._xp.where(has_image, n_non_nan_model_positions, 1) + ) ** len(self.data) + + chi_squared = -2.0 * ( -self._xp.log(n_permutations) + self._xp.sum(self.all_permutations_log_likelihoods()) ) + return self._xp.where(has_image, chi_squared, no_image_chi_squared) + class FitPositionsImagePairAllSolved(SolvedCentre, FitPositionsImagePairAll): """ diff --git a/autolens/point/solver/point_solver.py b/autolens/point/solver/point_solver.py index 8c28cf79e..86ab41eb0 100644 --- a/autolens/point/solver/point_solver.py +++ b/autolens/point/solver/point_solver.py @@ -123,36 +123,56 @@ def solve( ) # When no triangle traces to the source-plane coordinate -- e.g. the coordinate lies - # outside the region the image-plane grid tiles -- `filtered_means` is empty. The + # outside the region the image-plane grid tiles -- `filtered_means` comes back with shape + # (0, 2). The `[pair for pair in ...]` comprehension below then yields `[]`, and the # `Grid2DIrregular` built from it exposes `.array` as a bare list rather than a 2D array, - # so the `axis=1` reductions below used to raise + # so the `axis=1` reductions raised # `numpy.exceptions.AxisError: axis 1 is out of bounds for array of dimension 1`. # - # Zero images is a legitimate answer, not a failure, so return a correctly-shaped empty - # grid. `filtered_means` is a Python list, so its length is static (including under a - # `jax.jit` trace) and this branch is trace-safe. The JAX path pads to - # `MAX_CONTAINING_SIZE` and never reaches it. + # Zero images is a legitimate answer, not a failure, so build a correctly-shaped empty + # grid instead. + # + # `_filter_low_magnification` returns an `xp` array (its `List[Tuple[float, float]]` + # annotation is inaccurate), so `len(...)` is a leading-dimension query. Array shapes are + # static under a `jax.jit` trace, so this reads no traced *value* and the branch is + # trace-safe. In practice the JAX path builds its triangles at a padded fixed size + # (`MAX_CONTAINING_SIZE`), so it is the NumPy path that reaches this. if len(filtered_means) == 0: - logger.warning( - f"PointSolver.solve found no images for source-plane coordinate " - f"{tuple(source_plane_coordinate)}, so an empty grid is returned. This usually " - f"means the coordinate lies outside the region traced by the image-plane grid, or " - f"that every candidate image was removed by the magnification threshold." - ) + solution = xp.zeros((0, 2)) + + else: - return aa.Grid2DIrregular(xp.zeros((0, 2))) + solution = aa.Grid2DIrregular( + [pair for pair in filtered_means], xp=xp + ).array - solution = aa.Grid2DIrregular( - [pair for pair in filtered_means], xp=xp - ).array + is_nan = xp.isnan(solution).any(axis=1) + sentinel = xp.full_like(solution[0], fill_value=xp.inf) + solution = xp.where(is_nan[:, None], sentinel, solution) - is_nan = xp.isnan(solution).any(axis=1) - sentinel = xp.full_like(solution[0], fill_value=xp.inf) - solution = xp.where(is_nan[:, None], sentinel, solution) + if remove_infinities: - if remove_infinities: + solution = solution[~xp.isinf(solution).any(axis=1)] - solution = solution[~xp.isinf(solution).any(axis=1)] + # Warn on the *final* result rather than only on the branch above, because there are two + # distinct routes to an empty answer and only one goes through it: + # + # 1. no triangle contained the coordinate -> `filtered_means` is already length 0 + # 2. every candidate failed the magnification threshold -> `_filter_low_magnification` + # preserves the length and writes NaN rows, which become `inf` and are then stripped + # by `remove_infinities`, landing here at length 0 + # + # Only reachable on the NumPy path: the JAX path keeps its padded static shape, so + # `len(solution)` is non-zero there and this reads no traced value. + if len(solution) == 0: + + logger.warning( + f"PointSolver.solve found no images for source-plane coordinate " + f"{tuple(source_plane_coordinate)}, so an empty grid is returned. This means " + f"either that the coordinate lies outside the region traced by the image-plane " + f"grid, or that every candidate image was rejected by " + f"`magnification_threshold` (currently {self.magnification_threshold})." + ) return aa.Grid2DIrregular(solution) diff --git a/autolens/point/solver/shape_solver.py b/autolens/point/solver/shape_solver.py index 161db3bf4..08df41e80 100644 --- a/autolens/point/solver/shape_solver.py +++ b/autolens/point/solver/shape_solver.py @@ -304,7 +304,7 @@ def solve_triangles( would require {self.n_steps} subdivision steps. Decrease `pixel_scale_precision` (e.g. to {self.scale / 10.0:.3g} or smaller) so the - solver can resolve an image. + solver can refine the tiling. """ ) diff --git a/test_autolens/point/fit/positions/image/test_pair_all.py b/test_autolens/point/fit/positions/image/test_pair_all.py index bd4601169..4017f5b79 100644 --- a/test_autolens/point/fit/positions/image/test_pair_all.py +++ b/test_autolens/point/fit/positions/image/test_pair_all.py @@ -172,3 +172,60 @@ def test__fit_positions_image_pair_all_solved__source_plane_coordinate_feeds_sol # fixed regardless of the solved centre): matches the plain FitPositionsImagePairAll # value from the fixture-equivalent test above. assert fit.chi_squared == -2.0 * -4.40375330990644 + + +def test__no_model_positions__finite_no_image_floor_matching_siblings(data, noise_map): + """ + Regression: with no model positions, `n_permutations` is 0, so `-log(0)` is `+inf` while the + permutation sum is `-inf` and the two combined to NaN. + + `fitness.py` converts a NaN log-likelihood into `resample_figure_of_merit`, so the model was + silently resampled rather than scored -- exactly what the `no_image_residual` floor exists to + prevent. `FitPositionsImagePair` and `FitPositionsImagePairRepeat` already applied that floor; + `FitPositionsImagePairAll` did not. + + Surfaced by @rhayes777's PyAutoLens#531 once `PointSolver.solve` began returning an empty grid + instead of raising. + """ + + no_positions = al.Grid2DIrregular(np.zeros((0, 2))) + + kwargs = dict( + name="point_0", + data=data, + noise_map=noise_map, + tracer=tracer, + solver=al.mock.MockPointSolver(no_positions), + ) + + fit_all = al.FitPositionsImagePairAll(**kwargs) + + # `log_likelihood` is not asserted here: the module's `noise_map` fixture is a raw ndarray, + # which `fit_dataset.noise_normalization` cannot consume. That predates this fix and is why + # every test in this file asserts `chi_squared` only. `log_likelihood` finiteness is covered + # end-to-end against a real solver and an `ArrayIrregular` noise-map. + assert np.isfinite(fit_all.chi_squared) + + # the floor is on the same scale as the siblings', not merely finite + assert fit_all.chi_squared == pytest.approx( + al.FitPositionsImagePair(**kwargs).chi_squared + ) + assert fit_all.chi_squared == pytest.approx( + al.FitPositionsImagePairRepeat(**kwargs).chi_squared + ) + + +def test__model_positions_present__chi_squared_unchanged(data, noise_map): + """The no-image branch must not perturb the ordinary path.""" + + fit = al.FitPositionsImagePairAll( + name="point_0", + data=data, + noise_map=noise_map, + tracer=tracer, + solver=al.mock.MockPointSolver( + al.Grid2DIrregular([(-1.0749, -1.1), (1.19117, 1.175)]) + ), + ) + + assert fit.chi_squared == -2.0 * -4.40375330990644 diff --git a/test_autolens/point/triangles/test_solver_edge_cases.py b/test_autolens/point/triangles/test_solver_edge_cases.py index 0dc2f8fd7..46f3972f8 100644 --- a/test_autolens/point/triangles/test_solver_edge_cases.py +++ b/test_autolens/point/triangles/test_solver_edge_cases.py @@ -75,6 +75,53 @@ def test__source_outside_tiled_region__warns(solver_grid, lens_galaxy, caplog): assert "found no images" in caplog.text +def test__all_candidates_rejected_by_magnification__returns_empty_and_warns( + solver_grid, lens_galaxy, caplog +): + """ + The second route to an empty answer, which does NOT go through the length-0 branch. + + `_filter_low_magnification` preserves the array length and writes NaN rows, which become + `inf` and are stripped by `remove_infinities`. An earlier revision warned only on the + length-0 branch, so this route returned an empty grid silently. + """ + + source_plane_coordinate = (0.05, 0.02) + + solver = al.PointSolver.for_grid( + grid=solver_grid, + pixel_scale_precision=0.001, + magnification_threshold=1.0e100, + ) + + with caplog.at_level("WARNING"): + result = solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + assert len(result) == 0 + assert "found no images" in caplog.text + assert "magnification_threshold" in caplog.text + + +def test__successful_solve__does_not_warn(solver_grid, lens_galaxy, caplog): + """The warning must not fire on the ordinary path.""" + + source_plane_coordinate = (0.05, 0.02) + + solver = al.PointSolver.for_grid(grid=solver_grid, pixel_scale_precision=0.001) + + with caplog.at_level("WARNING"): + result = solver.solve( + tracer=_tracer(lens_galaxy, source_plane_coordinate), + source_plane_coordinate=source_plane_coordinate, + ) + + assert len(result) == 4 + assert "found no images" not in caplog.text + + @pytest.mark.parametrize("pixel_scale_precision", [0.05, 0.1, 0.2, 0.5]) def test__precision_coarser_than_triangle_scale__raises_naming_the_parameter( solver_grid, lens_galaxy, pixel_scale_precision