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 dfb3e393c..86ab41eb0 100644 --- a/autolens/point/solver/point_solver.py +++ b/autolens/point/solver/point_solver.py @@ -122,16 +122,57 @@ def solve( tracer=tracer, points=kept_triangles.means, xp=xp ) - 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) - - if remove_infinities: - - solution = solution[~xp.isinf(solution).any(axis=1)] + # When no triangle traces to the source-plane coordinate -- e.g. the coordinate lies + # 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 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 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: + + solution = xp.zeros((0, 2)) + + else: + + 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) + + if remove_infinities: + + 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 7b155b3dc..08df41e80 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 refine the tiling. + """ ) steps = list( 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 new file mode 100644 index 000000000..46f3972f8 --- /dev/null +++ b/test_autolens/point/triangles/test_solver_edge_cases.py @@ -0,0 +1,168 @@ +""" +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 + + +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 +): + """ + `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