fix: PointSolver zero-image and coarse-precision crashes, plus two consumer gaps

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)

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):
     """

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)

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(

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

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