diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py
index 3a95a1eee..5f2b2e523 100644
--- a/autolens/analysis/result.py
+++ b/autolens/analysis/result.py
@@ -31,7 +31,7 @@
 )
 from autolens.lens.tracer import Tracer
 from autolens.point.solver import PointSolver
-from autoconf.test_mode import skip_checks
+from autoconf.test_mode import is_test_mode, skip_checks
 
 logger = logging.getLogger(__name__)
 
@@ -301,6 +301,19 @@ def positions_likelihood_from(
             system is being analysed where the specific plane via its redshift is required to define whihch source
             galaxy is used to compute the multiple images.
 
+        Notes
+        -----
+        Test-mode safeguard (``PYAUTO_TEST_MODE``): integration tests intentionally fit
+        random / unphysical mass models and the resulting tracer often back-traces to
+        zero, NaN, or inf image-plane positions. Computing a position threshold from
+        such positions raises ``ValueError: zero-size array to reduction operation
+        fmax`` (or propagates NaN) inside ``positions_threshold_from``. When test mode
+        is active and the resolved positions are empty, contain NaN, or contain inf,
+        we substitute the synthetic pair ``[(1.0, 0.0), (-1.0, 0.0)]`` so the threshold
+        and ``PositionsLH`` still build cleanly and the script runs end-to-end.
+        Outside test mode this branch is never taken — bad positions still surface as
+        the original error, so production fits are not silently masked.
+
         Returns
         -------
         The `PositionsLH` object used to apply a likelihood penalty or resample the positions.
@@ -330,6 +343,15 @@ def positions_likelihood_from(
             np.asarray(positions.array if hasattr(positions, "array") else positions)
         )
 
+        if is_test_mode():
+            arr = positions.array
+            if arr.shape[0] < 2 or np.isnan(arr).any() or np.isinf(arr).any():
+                logger.warning(
+                    "positions_likelihood_from: empty/NaN/inf positions in PYAUTO_TEST_MODE — "
+                    "substituting synthetic fallback [(1.0, 0.0), (-1.0, 0.0)]."
+                )
+                positions = aa.Grid2DIrregular(values=[(1.0, 0.0), (-1.0, 0.0)])
+
         threshold = self.positions_threshold_from(
             factor=factor,
             minimum_threshold=minimum_threshold,
diff --git a/autolens/point/solver/point_solver.py b/autolens/point/solver/point_solver.py
index 29bb19f07..526823680 100644
--- a/autolens/point/solver/point_solver.py
+++ b/autolens/point/solver/point_solver.py
@@ -17,6 +17,7 @@
 default but can be retained for use inside a ``jax.jit``-traced function.
 """
 import logging
+import os
 from typing import Tuple, Optional
 
 import numpy as np
@@ -73,7 +74,22 @@ def solve(
         -------
         A ``Grid2DIrregular`` of image-plane coordinates, always numpy-backed even when the
         solver uses a JAX backend internally.
+
+        Notes
+        -----
+        Smoke-test short-circuit (``PYAUTO_SMALL_DATASETS``): the triangle-tiling solve
+        is the dominant cost in many simulator scripts and is meaningless on the
+        downsized grids used for fast smoke tests. When ``PYAUTO_SMALL_DATASETS=1`` is
+        set the solver returns the fixed pair ``[(1.0, 0.0), (0.0, 1.0)]`` immediately,
+        skipping ``solve_triangles`` entirely. The two coordinates are well separated
+        so any downstream ``positions_likelihood_from`` / threshold calculation behaves
+        normally. ``PYAUTO_SMALL_DATASETS`` is a smoke-test-only flag and is never set
+        inside a ``jax.jit`` trace, so a plain numpy-backed ``Grid2DIrregular`` is safe
+        here even when the surrounding analysis uses ``xp=jnp``.
         """
+        if os.environ.get("PYAUTO_SMALL_DATASETS") == "1":
+            return aa.Grid2DIrregular(values=[(1.0, 0.0), (0.0, 1.0)])
+
         kept_triangles = super().solve_triangles(
             tracer=tracer,
             shape=Point(*source_plane_coordinate),
diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py
index a5835ca1c..4d99fd7d5 100644
--- a/test_autolens/analysis/test_result.py
+++ b/test_autolens/analysis/test_result.py
@@ -248,6 +248,38 @@ def test__positions_likelihood_from(analysis_imaging_7x7):
     assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
 
 
+def test__positions_likelihood_from__test_mode_fallback(
+    monkeypatch, analysis_imaging_7x7,
+):
+    monkeypatch.setenv("PYAUTO_TEST_MODE", "2")
+
+    tracer = al.Tracer(
+        galaxies=[
+            al.Galaxy(
+                redshift=0.5,
+                mass=al.mp.Isothermal(
+                    centre=(0.1, 0.0), einstein_radius=1.0, ell_comps=(0.0, 0.0)
+                ),
+            ),
+            al.Galaxy(redshift=1.0, bulge=al.lp.SersicSph(centre=(0.0, 0.0))),
+        ]
+    )
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    empty_positions = al.Grid2DIrregular(np.empty((0, 2)))
+
+    positions_likelihood = result.positions_likelihood_from(
+        factor=0.1, minimum_threshold=0.2, positions=empty_positions
+    )
+
+    assert isinstance(positions_likelihood, al.PositionsLH)
+    assert len(positions_likelihood.positions) == 2
+    assert positions_likelihood.positions[0] == pytest.approx((1.0, 0.0))
+    assert positions_likelihood.positions[1] == pytest.approx((-1.0, 0.0))
+
+
 def test__positions_likelihood_from__mass_centre_radial_distance_min(
     analysis_imaging_7x7,
 ):