diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py
index 3e65aac16..9dcfcd3f6 100644
--- a/autolens/analysis/result.py
+++ b/autolens/analysis/result.py
@@ -1,527 +1,523 @@
-"""
-Result classes for **PyAutoLens** model fits.
-
-This module defines the ``Result`` and ``ResultDataset`` classes that wrap the output of
-a ``PyAutoFit`` non-linear search for lensing model fits.  Key properties include:
-
-- ``max_log_likelihood_tracer`` — reconstructs the best-fit ``Tracer`` from the stored
-  model instance.
-- ``source_plane_light_profile_centre_from`` — extracts the source-plane light-profile
-  centre, used by automatic position updating to seed position priors for the next
-  search in a SLaM pipeline.
-- ``image_plane_multiple_image_positions_from`` — uses a ``PointSolver`` to compute the
-  image-plane multiple-image positions predicted by the best-fit tracer, optionally
-  filtered by a minimum flux threshold.
-
-These results feed directly into downstream pipeline stages and post-processing scripts.
-"""
-import json
-import logging
-import os
-import numpy as np
-from typing import List, Optional, Union
-
-import autoarray as aa
-import autogalaxy as ag
-
-from autogalaxy.analysis.result import ResultDataset as AgResultDataset
-
-from autolens.analysis.positions import PositionsLH
-from autolens.point.max_separation import (
-    SourceMaxSeparation,
-)
-from autolens.lens.tracer import Tracer
-from autolens.point.solver import PointSolver
-from autoconf.test_mode import is_test_mode, skip_checks
-
-logger = logging.getLogger(__name__)
-
-
-class Result(AgResultDataset):
-    @property
-    def max_log_likelihood_tracer(self) -> Tracer:
-        """
-        An instance of a `Tracer` corresponding to the maximum log likelihood model inferred by the non-linear search.
-        """
-        return self.analysis.tracer_via_instance_from(instance=self.instance)
-
-    def source_plane_light_profile_centre_from(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        Return a light profile centre of one of the a galaxies in the maximum log likelihood `Tracer`'s source-plane.
-        If there are multiple light profiles, the first light profile's centre is returned.
-
-        These centres are used by automatic position updating to determine the best-fit lens model's image-plane
-        multiple-image positions.
-        """
-        try:
-            plane_index = self.max_log_likelihood_tracer.plane_index_via_redshift_from(
-                redshift=plane_redshift
-            )
-        except TypeError:
-            plane_index = -1
-
-        centre = self.max_log_likelihood_tracer.planes[plane_index].extract_attribute(
-            cls=ag.LightProfile, attr_name="centre"
-        )
-        if centre is not None:
-            return aa.Grid2DIrregular(values=[np.asarray(centre[0])])
-
-    def source_plane_centre_from(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        Return the centre of a source-plane galaxy via the following criteria:
-
-        1) If the source plane contains only light profiles, return the first light's centre.
-        2) If it contains an `LEq` return the centre of its brightest pixel instead.
-
-        These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model
-        (and thus tracer) by back-tracing the centres to the image plane via the mass model.
-        """
-        return self.source_plane_light_profile_centre_from(
-            plane_redshift=plane_redshift
-        )
-
-    def image_plane_multiple_image_positions(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        Backwards ray-trace the source-plane centres (see above) to the image-plane via the mass model, to determine
-        the multiple image position of the source(s) in the image-plane.
-
-        These image-plane positions are used by the next search in a pipeline if automatic position updating is turned
-        on.
-
-        Parameters
-        ----------
-        plane_redshift
-            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
-            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.
-        """
-
-        grid = self.analysis.dataset.mask.derive_grid.all_false
-
-        solver = PointSolver.for_grid(
-            grid=grid, pixel_scale_precision=0.001
-        )
-
-        source_plane_centre = self.source_plane_centre_from(
-            plane_redshift=plane_redshift
-        )
-
-        multiple_images = solver.solve(
-            tracer=self.max_log_likelihood_tracer,
-            source_plane_coordinate=source_plane_centre.in_list[0],
-            xp=self.analysis._xp,
-            plane_redshift=plane_redshift,
-        )
-
-        if multiple_images.shape[0] == 1:
-            return self.image_plane_multiple_image_positions_for_single_image_from()
-
-        return aa.Grid2DIrregular(values=multiple_images)
-
-    def image_plane_multiple_image_positions_for_single_image_from(
-        self, plane_redshift: Optional[float] = None, increments: int = 20
-    ) -> aa.Grid2DIrregular:
-        """
-        If the standard point solver only locates one multiple image, finds one or more additional images, which are
-        not technically multiple image in the point source regime, but are close enough to it they can be used
-        in a position threshold likelihood.
-
-        This is performed by incrementally moving the source-plane centre's coordinates towards the centre of the
-        source-plane at (0.0", 0.0"). This ensures that the centre will eventually go inside the caustic, where
-        multiple images are formed.
-
-        To move the source-plane centre, the original source-plane centre is multiplied by a factor that decreases
-        from 1.0 to 0.0 in increments of 1/increments. For example, if the source-plane centre is (1.0", -0.5") and
-        the `factor` is 0.5, the input source-plane centre is (0.5", -0.25").
-
-        The multiple images are always computed for the same mass model, thus they will always be valid multiple images
-        for the model being fitted, but as the factor decrease the multiple images may move furhter from their observed
-        positions.
-
-        Parameters
-        ----------
-        plane_redshift
-            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
-            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.
-        increments
-            The number of increments the source-plane centre is moved to compute multiple images.
-        """
-
-        logger.info(
-            """
-        Could not find multiple images for maximum likelihood lens model.
-        
-        Incrementally moving source centre inwards towards centre of source-plane until caustic crossing occurs
-        and multiple images are formed.        
-        """
-        )
-
-        grid = self.analysis.dataset.mask.derive_grid.all_false
-
-        centre = self.source_plane_centre_from(plane_redshift=plane_redshift).in_list[0]
-
-        solver = PointSolver.for_grid(
-            grid=grid, pixel_scale_precision=0.001
-        )
-
-        for i in range(1, increments):
-            factor = 1.0 - (1.0 * (i / increments))
-
-            multiple_images = solver.solve(
-                tracer=self.max_log_likelihood_tracer,
-                source_plane_coordinate=(centre[0] * factor, centre[1] * factor),
-                xp=self.analysis._xp,
-                plane_redshift=plane_redshift,
-            )
-
-            if multiple_images.shape[0] > 1:
-                return aa.Grid2DIrregular(values=multiple_images)
-
-        logger.info(
-            """
-        Could not find multiple images for maximum likelihood lens model, even after incrementally moving the source
-        centre inwards to the centre of the source-plane.
-
-        Set the multiple image postiions to two images at (1.0", 1.0") so code continues to run.
-        """
-        )
-        return aa.Grid2DIrregular(values=[(1.0, 1.0), (1.0, 1.0)])
-
-    def positions_threshold_from(
-        self,
-        factor=1.0,
-        minimum_threshold=None,
-        maximum_threshold=None,
-        positions: Optional[aa.Grid2DIrregular] = None,
-        plane_redshift: Optional[float] = None,
-    ) -> float:
-        """
-        Compute a new position threshold from these results corresponding to the image-plane multiple image positions
-        of the maximum log likelihood `Tracer` ray-traced to the source-plane.
-
-        First, we ray-trace forward the multiple-image's to the source-plane via the mass model to determine how far
-        apart they are separated. We take the maximum source-plane separation of these points and multiple this by
-        the auto_positions_factor to determine a new positions threshold. This value may also be rounded up to the
-        input `auto_positions_minimum_threshold`.
-
-        This is used for non-linear search chaining, specifically updating the position threshold of a new model-fit
-        using the maximum likelihood model of a previous search.
-
-        The above behaviour assumes a single lens and single source plane, if there are multiple source planes (e.g.
-        a double source plane lens) the input `plane_redshift_list` can be used to specify which source planea
-        the position threshold is computed for.
-
-        Parameters
-        ----------
-        factor
-            The value the computed threshold is multiplied by to make the position threshold larger or smaller than the
-            maximum log likelihood model's threshold.
-        minimum_threshold
-            The output threshold is rounded up to this value if it is below it, to avoid extremely small threshold
-            values.
-        maximum_threshold
-            The output threshold is rounded down to this value if it is above it, to avoid extremely large threshold
-            values.
-        positions
-            If input, these positions are used instead of the computed multiple image positions from the lens mass
-            model.
-        plane_redshift
-            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
-            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.
-
-        Returns
-        -------
-        float
-            The maximum source plane separation of this results maximum likelihood `Tracer` multiple images multiplied
-            by `factor` and rounded up to the `threshold`.
-        """
-
-        tracer = Tracer(galaxies=self.max_log_likelihood_galaxies)
-
-        positions = (
-            self.image_plane_multiple_image_positions(plane_redshift=plane_redshift)
-            if positions is None
-            else positions
-        )
-
-        positions_fits = SourceMaxSeparation(
-            data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift
-        )
-
-        threshold = factor * np.nanmax(positions_fits.max_separation_of_plane_positions)
-
-        if minimum_threshold is not None:
-            if threshold < minimum_threshold:
-                return minimum_threshold
-
-        if maximum_threshold is not None:
-            if threshold > maximum_threshold:
-                return maximum_threshold
-
-        return threshold
-
-    def _cached_multiple_image_positions_from(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        The multiple image positions of the maximum log likelihood lens model, loaded
-        from this result's own ``files/`` folder when previously solved, else solved
-        via the point solver and persisted there.
-
-        Solving the multiple image positions runs a point-solver grid search over the
-        maximum log likelihood tracer, which on a resumed pipeline (e.g. SLaM) pays a
-        fresh JIT compile and dominates resume overhead (autolens_profiling#70) — the
-        solved positions themselves are a pure product of the completed fit, so they
-        are cached as ``files/multiple_image_positions[_plane_<z>].json``. Staleness
-        is structurally guarded: a changed model or search produces a new search
-        identifier and a fresh output directory with no cache file. Results with no
-        on-disk output (e.g. ``NullPaths``) always solve.
-        """
-        from pathlib import Path
-
-        name = "multiple_image_positions"
-        if plane_redshift is not None:
-            name += f"_plane_{str(plane_redshift).replace('.', '_')}"
-
-        files_path = getattr(getattr(self, "paths", None), "_files_path", None)
-        cache_path = (
-            Path(files_path) / f"{name}.json"
-            if files_path is not None and Path(files_path).is_dir()
-            else None
-        )
-
-        if cache_path is not None and cache_path.exists():
-            with open(cache_path) as f:
-                return aa.Grid2DIrregular(values=[tuple(p) for p in json.load(f)])
-
-        positions = self.image_plane_multiple_image_positions(
-            plane_redshift=plane_redshift
-        )
-
-        if cache_path is not None:
-            with open(cache_path, "w") as f:
-                json.dump(np.asarray(positions.array).tolist(), f)
-
-            # Preserve the cache in the search's zip — a resumed search's
-            # paths.restore() wipes the output dir and re-extracts the zip,
-            # destroying any file written only to files/ after completion.
-            from autogalaxy.analysis.adapt_images.adapt_images import (
-                _append_to_search_zip,
-            )
-
-            _append_to_search_zip(self.paths, cache_path)
-
-        return positions
-
-    def positions_likelihood_from(
-        self,
-        factor=1.0,
-        minimum_threshold=None,
-        maximum_threshold=None,
-        positions: Optional[aa.Grid2DIrregular] = None,
-        mass_centre_radial_distance_min: float = None,
-        plane_redshift: Optional[float] = None,
-    ) -> PositionsLH:
-        """
-        Returns a `PositionsLH` object from the result of a lens model-fit, where the maximum log likelihood mass
-        and source models are used to determine the multiple image positions in the image-plane and source-plane
-        and ray-trace them to the source-plane to determine the threshold.
-
-        In chained fits, for example the SLaM pipelines, this means that a simple initial fit (e.g. SIE mass model,
-        parametric source) can be used to determine the multiple image positions and threshold for a more complex
-        subsequent fit (e.g. power-law mass model, pixelized source).
-
-        The mass model central image is removed from the solution, as this is rarely physically observed and therefore
-        should not be included in the likelihood penalty or resampling. It is removed by setting a positive
-        magnification threshold in the `PointSolver`. For strange lens models the central image may still be
-        solved for, in which case the `mass_centre_radial_distance_min` parameter can be used to remove it.
-
-        Parameters
-        ----------
-        factor
-            The value the computed threshold is multiplied by to make the position threshold larger or smaller than the
-            maximum log likelihood model's threshold.
-        minimum_threshold
-            The output threshold is rounded up to this value if it is below it, to avoid extremely small threshold
-            values.
-        maximum_threshold
-            The output threshold is rounded down to this value if it is above it, to avoid extremely large threshold
-            values.
-        positions
-            If input, these positions are used instead of the computed multiple image positions from the lens mass
-            model.
-        mass_centre_radial_distance_min
-            The minimum radial distance from the mass model centre that a multiple image position must be to be
-            included in the likelihood penalty or resampling. If `None` all positions are used. This is an additional
-            method to remove central images that may make it through the point solver's magnification threshold.
-        plane_redshift
-            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
-            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.
-
-        Skip-checks safeguard (``PYAUTO_SKIP_CHECKS``): when validation checks are
-        skipped *and* test mode is active, we return a ``PositionsLH`` built from the
-        same synthetic ``[(1.0, 0.0), (-1.0, 0.0)]`` pair (with threshold set to
-        ``minimum_threshold`` or 0.5 if unset). This avoids the expensive
-        multiple-image-position computation while still giving callers a usable
-        ``PositionsLH``, so workspace scripts that chain ``.positions`` continue to
-        work end-to-end. With skip-checks alone (no test mode) the function still
-        returns ``None`` — the production opt-out behaviour is unchanged.
-
-        Returns
-        -------
-        The `PositionsLH` object used to apply a likelihood penalty or resample the positions.
-        """
-
-        if skip_checks():
-            if is_test_mode():
-                synthetic_positions = aa.Grid2DIrregular(
-                    values=[(1.0, 0.0), (-1.0, 0.0)]
-                )
-                synthetic_threshold = (
-                    minimum_threshold if minimum_threshold is not None else 0.5
-                )
-                return PositionsLH(
-                    positions=synthetic_positions,
-                    threshold=synthetic_threshold,
-                    plane_redshift=plane_redshift,
-                )
-            return
-
-        if positions is None:
-            positions = self._cached_multiple_image_positions_from(
-                plane_redshift=plane_redshift
-            )
-
-        if mass_centre_radial_distance_min is not None:
-            mass_centre = self.max_log_likelihood_tracer.extract_attribute(
-                cls=ag.mp.MassProfile, attr_name="centre"
-            )
-
-            distances = positions.distances_to_coordinate_from(
-                coordinate=mass_centre[0]
-            )
-
-            positions = positions[distances > mass_centre_radial_distance_min]
-
-        positions = aa.Grid2DIrregular(
-            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,
-            maximum_threshold=maximum_threshold,
-            positions=positions,
-            plane_redshift=plane_redshift,
-        )
-
-        if maximum_threshold is not None:
-            threshold = min(threshold, maximum_threshold)
-
-        return PositionsLH(
-            positions=positions, threshold=threshold, plane_redshift=plane_redshift
-        )
-
-
-class ResultDataset(Result):
-    @property
-    def max_log_likelihood_tracer(self) -> Tracer:
-        """
-        An instance of a `Tracer` corresponding to the maximum log likelihood model inferred by the non-linear search.
-
-        If a dataset is fitted the adapt images of the adapt image must first be associated with each galaxy.
-        """
-        instance = self.analysis.instance_with_associated_adapt_images_from(
-            instance=self.instance
-        )
-
-        return self.analysis.tracer_via_instance_from(instance=instance)
-
-    def source_plane_centre_from(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        Return the centre of a source-plane galaxy via the following criteria:
-
-        1) If the source plane contains only light profiles, return the first light's centre.
-        2) If it contains an `LEq` return the centre of its brightest pixel instead.
-
-        These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model
-        (and thus tracer) by back-tracing the centres to the image plane via the mass model.
-        """
-        source_plane_inversion_centre = self.source_plane_inversion_centre_from(
-            plane_redshift=plane_redshift
-        )
-
-        if source_plane_inversion_centre is not None:
-            return source_plane_inversion_centre
-
-        return self.source_plane_light_profile_centre_from(
-            plane_redshift=plane_redshift
-        )
-
-    def source_plane_inversion_centre_from(
-        self, plane_redshift: Optional[float] = None
-    ) -> aa.Grid2DIrregular:
-        """
-        Returns the centre of the brightest source pixel(s) of an `LEq`.
-
-        These centres are used by automatic position updating to determine the best-fit lens model's image-plane
-        multiple-image positions.
-        """
-
-        if self.max_log_likelihood_fit.inversion is not None:
-            if self.max_log_likelihood_fit.inversion.has(cls=aa.Mapper):
-
-                try:
-                    plane_index = (
-                        self.max_log_likelihood_tracer.plane_index_via_redshift_from(
-                            redshift=plane_redshift
-                        )
-                    )
-                    # plane_indexes_with_pixelizations is a list of plane indices that
-                    # have a pixelization, in mapper order. The mapper index for a given
-                    # plane is the position of that plane index within the list, not the
-                    # element at [plane_index].
-                    mapper_index = (
-                        self.max_log_likelihood_tracer.plane_indexes_with_pixelizations.index(
-                            plane_index
-                        )
-                    )
-                except (TypeError, ValueError):
-                    mapper_index = 0
-
-                inversion = self.max_log_likelihood_fit.inversion
-
-                return inversion.max_pixel_centre(mapper_index=mapper_index)
+"""
+Result classes for **PyAutoLens** model fits.
+
+This module defines the ``Result`` and ``ResultDataset`` classes that wrap the output of
+a ``PyAutoFit`` non-linear search for lensing model fits.  Key properties include:
+
+- ``max_log_likelihood_tracer`` — reconstructs the best-fit ``Tracer`` from the stored
+  model instance.
+- ``source_plane_light_profile_centre_from`` — extracts the source-plane light-profile
+  centre, used by automatic position updating to seed position priors for the next
+  search in a SLaM pipeline.
+- ``image_plane_multiple_image_positions_from`` — uses a ``PointSolver`` to compute the
+  image-plane multiple-image positions predicted by the best-fit tracer, optionally
+  filtered by a minimum flux threshold.
+
+These results feed directly into downstream pipeline stages and post-processing scripts.
+"""
+import json
+import logging
+import os
+import numpy as np
+from typing import List, Optional, Union
+
+import autoarray as aa
+import autogalaxy as ag
+
+from autogalaxy.analysis.result import ResultDataset as AgResultDataset
+
+from autolens.analysis.positions import PositionsLH
+from autolens.point.max_separation import (
+    SourceMaxSeparation,
+)
+from autolens.lens.tracer import Tracer
+from autolens.point.solver import PointSolver
+from autoconf.test_mode import is_test_mode, skip_checks
+
+logger = logging.getLogger(__name__)
+
+
+class Result(AgResultDataset):
+    @property
+    def max_log_likelihood_tracer(self) -> Tracer:
+        """
+        An instance of a `Tracer` corresponding to the maximum log likelihood model inferred by the non-linear search.
+        """
+        return self.analysis.tracer_via_instance_from(instance=self.instance)
+
+    def source_plane_light_profile_centre_from(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        Return a light profile centre of one of the a galaxies in the maximum log likelihood `Tracer`'s source-plane.
+        If there are multiple light profiles, the first light profile's centre is returned.
+
+        These centres are used by automatic position updating to determine the best-fit lens model's image-plane
+        multiple-image positions.
+        """
+        try:
+            plane_index = self.max_log_likelihood_tracer.plane_index_via_redshift_from(
+                redshift=plane_redshift
+            )
+        except TypeError:
+            plane_index = -1
+
+        centre = self.max_log_likelihood_tracer.planes[plane_index].extract_attribute(
+            cls=ag.LightProfile, attr_name="centre"
+        )
+        if centre is not None:
+            return aa.Grid2DIrregular(values=[np.asarray(centre[0])])
+
+    def source_plane_centre_from(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        Return the centre of a source-plane galaxy via the following criteria:
+
+        1) If the source plane contains only light profiles, return the first light's centre.
+        2) If it contains an `LEq` return the centre of its brightest pixel instead.
+
+        These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model
+        (and thus tracer) by back-tracing the centres to the image plane via the mass model.
+        """
+        return self.source_plane_light_profile_centre_from(
+            plane_redshift=plane_redshift
+        )
+
+    def image_plane_multiple_image_positions(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        Backwards ray-trace the source-plane centres (see above) to the image-plane via the mass model, to determine
+        the multiple image position of the source(s) in the image-plane.
+
+        These image-plane positions are used by the next search in a pipeline if automatic position updating is turned
+        on.
+
+        Parameters
+        ----------
+        plane_redshift
+            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
+            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.
+        """
+
+        grid = self.analysis.dataset.mask.derive_grid.all_false
+
+        solver = PointSolver.for_grid(
+            grid=grid, pixel_scale_precision=0.001
+        )
+
+        source_plane_centre = self.source_plane_centre_from(
+            plane_redshift=plane_redshift
+        )
+
+        multiple_images = solver.solve(
+            tracer=self.max_log_likelihood_tracer,
+            source_plane_coordinate=source_plane_centre.in_list[0],
+            xp=self.analysis._xp,
+            plane_redshift=plane_redshift,
+        )
+
+        if multiple_images.shape[0] == 1:
+            return self.image_plane_multiple_image_positions_for_single_image_from()
+
+        return aa.Grid2DIrregular(values=multiple_images)
+
+    def image_plane_multiple_image_positions_for_single_image_from(
+        self, plane_redshift: Optional[float] = None, increments: int = 20
+    ) -> aa.Grid2DIrregular:
+        """
+        If the standard point solver only locates one multiple image, finds one or more additional images, which are
+        not technically multiple image in the point source regime, but are close enough to it they can be used
+        in a position threshold likelihood.
+
+        This is performed by incrementally moving the source-plane centre's coordinates towards the centre of the
+        source-plane at (0.0", 0.0"). This ensures that the centre will eventually go inside the caustic, where
+        multiple images are formed.
+
+        To move the source-plane centre, the original source-plane centre is multiplied by a factor that decreases
+        from 1.0 to 0.0 in increments of 1/increments. For example, if the source-plane centre is (1.0", -0.5") and
+        the `factor` is 0.5, the input source-plane centre is (0.5", -0.25").
+
+        The multiple images are always computed for the same mass model, thus they will always be valid multiple images
+        for the model being fitted, but as the factor decrease the multiple images may move furhter from their observed
+        positions.
+
+        Parameters
+        ----------
+        plane_redshift
+            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
+            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.
+        increments
+            The number of increments the source-plane centre is moved to compute multiple images.
+        """
+
+        logger.info(
+            """
+        Could not find multiple images for maximum likelihood lens model.
+        
+        Incrementally moving source centre inwards towards centre of source-plane until caustic crossing occurs
+        and multiple images are formed.        
+        """
+        )
+
+        grid = self.analysis.dataset.mask.derive_grid.all_false
+
+        centre = self.source_plane_centre_from(plane_redshift=plane_redshift).in_list[0]
+
+        solver = PointSolver.for_grid(
+            grid=grid, pixel_scale_precision=0.001
+        )
+
+        for i in range(1, increments):
+            factor = 1.0 - (1.0 * (i / increments))
+
+            multiple_images = solver.solve(
+                tracer=self.max_log_likelihood_tracer,
+                source_plane_coordinate=(centre[0] * factor, centre[1] * factor),
+                xp=self.analysis._xp,
+                plane_redshift=plane_redshift,
+            )
+
+            if multiple_images.shape[0] > 1:
+                return aa.Grid2DIrregular(values=multiple_images)
+
+        logger.info(
+            """
+        Could not find multiple images for maximum likelihood lens model, even after incrementally moving the source
+        centre inwards to the centre of the source-plane.
+
+        Set the multiple image postiions to two images at (1.0", 1.0") so code continues to run.
+        """
+        )
+        return aa.Grid2DIrregular(values=[(1.0, 1.0), (1.0, 1.0)])
+
+    def positions_threshold_from(
+        self,
+        factor=1.0,
+        minimum_threshold=None,
+        maximum_threshold=None,
+        positions: Optional[aa.Grid2DIrregular] = None,
+        plane_redshift: Optional[float] = None,
+    ) -> float:
+        """
+        Compute a new position threshold from these results corresponding to the image-plane multiple image positions
+        of the maximum log likelihood `Tracer` ray-traced to the source-plane.
+
+        First, we ray-trace forward the multiple-image's to the source-plane via the mass model to determine how far
+        apart they are separated. We take the maximum source-plane separation of these points and multiple this by
+        the auto_positions_factor to determine a new positions threshold. This value may also be rounded up to the
+        input `auto_positions_minimum_threshold`.
+
+        This is used for non-linear search chaining, specifically updating the position threshold of a new model-fit
+        using the maximum likelihood model of a previous search.
+
+        The above behaviour assumes a single lens and single source plane, if there are multiple source planes (e.g.
+        a double source plane lens) the input `plane_redshift_list` can be used to specify which source planea
+        the position threshold is computed for.
+
+        Parameters
+        ----------
+        factor
+            The value the computed threshold is multiplied by to make the position threshold larger or smaller than the
+            maximum log likelihood model's threshold.
+        minimum_threshold
+            The output threshold is rounded up to this value if it is below it, to avoid extremely small threshold
+            values.
+        maximum_threshold
+            The output threshold is rounded down to this value if it is above it, to avoid extremely large threshold
+            values.
+        positions
+            If input, these positions are used instead of the computed multiple image positions from the lens mass
+            model.
+        plane_redshift
+            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
+            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.
+
+        Returns
+        -------
+        float
+            The maximum source plane separation of this results maximum likelihood `Tracer` multiple images multiplied
+            by `factor` and rounded up to the `threshold`.
+        """
+
+        tracer = Tracer(galaxies=self.max_log_likelihood_galaxies)
+
+        positions = (
+            self.image_plane_multiple_image_positions(plane_redshift=plane_redshift)
+            if positions is None
+            else positions
+        )
+
+        positions_fits = SourceMaxSeparation(
+            data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift
+        )
+
+        threshold = factor * np.nanmax(positions_fits.max_separation_of_plane_positions)
+
+        if minimum_threshold is not None:
+            if threshold < minimum_threshold:
+                return minimum_threshold
+
+        if maximum_threshold is not None:
+            if threshold > maximum_threshold:
+                return maximum_threshold
+
+        return threshold
+
+    def _cached_multiple_image_positions_from(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        The multiple image positions of the maximum log likelihood lens model, loaded
+        from this result's own ``files/`` folder when previously solved, else solved
+        via the point solver and persisted there.
+
+        Solving the multiple image positions runs a point-solver grid search over the
+        maximum log likelihood tracer, which on a resumed pipeline (e.g. SLaM) pays a
+        fresh JIT compile and dominates resume overhead (autolens_profiling#70) — the
+        solved positions themselves are a pure product of the completed fit, so they
+        are cached as ``files/multiple_image_positions[_plane_<z>].json``. Staleness
+        is structurally guarded: a changed model or search produces a new search
+        identifier and a fresh output directory with no cache file. Results with no
+        on-disk output (e.g. ``NullPaths``) always solve.
+        """
+        from pathlib import Path
+
+        name = "multiple_image_positions"
+        if plane_redshift is not None:
+            name += f"_plane_{str(plane_redshift).replace('.', '_')}"
+
+        files_path = getattr(getattr(self, "paths", None), "_files_path", None)
+        cache_path = (
+            Path(files_path) / f"{name}.json"
+            if files_path is not None and Path(files_path).is_dir()
+            else None
+        )
+
+        if cache_path is not None and cache_path.exists():
+            with open(cache_path) as f:
+                return aa.Grid2DIrregular(values=[tuple(p) for p in json.load(f)])
+
+        positions = self.image_plane_multiple_image_positions(
+            plane_redshift=plane_redshift
+        )
+
+        if cache_path is not None:
+            with open(cache_path, "w") as f:
+                json.dump(np.asarray(positions.array).tolist(), f)
+
+            # Preserve the cache in the search's zip — a resumed search's
+            # paths.restore() wipes the output dir and re-extracts the zip,
+            # destroying any file written only to files/ after completion.
+            self.paths.preserve_in_zip(cache_path)
+
+        return positions
+
+    def positions_likelihood_from(
+        self,
+        factor=1.0,
+        minimum_threshold=None,
+        maximum_threshold=None,
+        positions: Optional[aa.Grid2DIrregular] = None,
+        mass_centre_radial_distance_min: float = None,
+        plane_redshift: Optional[float] = None,
+    ) -> PositionsLH:
+        """
+        Returns a `PositionsLH` object from the result of a lens model-fit, where the maximum log likelihood mass
+        and source models are used to determine the multiple image positions in the image-plane and source-plane
+        and ray-trace them to the source-plane to determine the threshold.
+
+        In chained fits, for example the SLaM pipelines, this means that a simple initial fit (e.g. SIE mass model,
+        parametric source) can be used to determine the multiple image positions and threshold for a more complex
+        subsequent fit (e.g. power-law mass model, pixelized source).
+
+        The mass model central image is removed from the solution, as this is rarely physically observed and therefore
+        should not be included in the likelihood penalty or resampling. It is removed by setting a positive
+        magnification threshold in the `PointSolver`. For strange lens models the central image may still be
+        solved for, in which case the `mass_centre_radial_distance_min` parameter can be used to remove it.
+
+        Parameters
+        ----------
+        factor
+            The value the computed threshold is multiplied by to make the position threshold larger or smaller than the
+            maximum log likelihood model's threshold.
+        minimum_threshold
+            The output threshold is rounded up to this value if it is below it, to avoid extremely small threshold
+            values.
+        maximum_threshold
+            The output threshold is rounded down to this value if it is above it, to avoid extremely large threshold
+            values.
+        positions
+            If input, these positions are used instead of the computed multiple image positions from the lens mass
+            model.
+        mass_centre_radial_distance_min
+            The minimum radial distance from the mass model centre that a multiple image position must be to be
+            included in the likelihood penalty or resampling. If `None` all positions are used. This is an additional
+            method to remove central images that may make it through the point solver's magnification threshold.
+        plane_redshift
+            The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens
+            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.
+
+        Skip-checks safeguard (``PYAUTO_SKIP_CHECKS``): when validation checks are
+        skipped *and* test mode is active, we return a ``PositionsLH`` built from the
+        same synthetic ``[(1.0, 0.0), (-1.0, 0.0)]`` pair (with threshold set to
+        ``minimum_threshold`` or 0.5 if unset). This avoids the expensive
+        multiple-image-position computation while still giving callers a usable
+        ``PositionsLH``, so workspace scripts that chain ``.positions`` continue to
+        work end-to-end. With skip-checks alone (no test mode) the function still
+        returns ``None`` — the production opt-out behaviour is unchanged.
+
+        Returns
+        -------
+        The `PositionsLH` object used to apply a likelihood penalty or resample the positions.
+        """
+
+        if skip_checks():
+            if is_test_mode():
+                synthetic_positions = aa.Grid2DIrregular(
+                    values=[(1.0, 0.0), (-1.0, 0.0)]
+                )
+                synthetic_threshold = (
+                    minimum_threshold if minimum_threshold is not None else 0.5
+                )
+                return PositionsLH(
+                    positions=synthetic_positions,
+                    threshold=synthetic_threshold,
+                    plane_redshift=plane_redshift,
+                )
+            return
+
+        if positions is None:
+            positions = self._cached_multiple_image_positions_from(
+                plane_redshift=plane_redshift
+            )
+
+        if mass_centre_radial_distance_min is not None:
+            mass_centre = self.max_log_likelihood_tracer.extract_attribute(
+                cls=ag.mp.MassProfile, attr_name="centre"
+            )
+
+            distances = positions.distances_to_coordinate_from(
+                coordinate=mass_centre[0]
+            )
+
+            positions = positions[distances > mass_centre_radial_distance_min]
+
+        positions = aa.Grid2DIrregular(
+            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,
+            maximum_threshold=maximum_threshold,
+            positions=positions,
+            plane_redshift=plane_redshift,
+        )
+
+        if maximum_threshold is not None:
+            threshold = min(threshold, maximum_threshold)
+
+        return PositionsLH(
+            positions=positions, threshold=threshold, plane_redshift=plane_redshift
+        )
+
+
+class ResultDataset(Result):
+    @property
+    def max_log_likelihood_tracer(self) -> Tracer:
+        """
+        An instance of a `Tracer` corresponding to the maximum log likelihood model inferred by the non-linear search.
+
+        If a dataset is fitted the adapt images of the adapt image must first be associated with each galaxy.
+        """
+        instance = self.analysis.instance_with_associated_adapt_images_from(
+            instance=self.instance
+        )
+
+        return self.analysis.tracer_via_instance_from(instance=instance)
+
+    def source_plane_centre_from(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        Return the centre of a source-plane galaxy via the following criteria:
+
+        1) If the source plane contains only light profiles, return the first light's centre.
+        2) If it contains an `LEq` return the centre of its brightest pixel instead.
+
+        These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model
+        (and thus tracer) by back-tracing the centres to the image plane via the mass model.
+        """
+        source_plane_inversion_centre = self.source_plane_inversion_centre_from(
+            plane_redshift=plane_redshift
+        )
+
+        if source_plane_inversion_centre is not None:
+            return source_plane_inversion_centre
+
+        return self.source_plane_light_profile_centre_from(
+            plane_redshift=plane_redshift
+        )
+
+    def source_plane_inversion_centre_from(
+        self, plane_redshift: Optional[float] = None
+    ) -> aa.Grid2DIrregular:
+        """
+        Returns the centre of the brightest source pixel(s) of an `LEq`.
+
+        These centres are used by automatic position updating to determine the best-fit lens model's image-plane
+        multiple-image positions.
+        """
+
+        if self.max_log_likelihood_fit.inversion is not None:
+            if self.max_log_likelihood_fit.inversion.has(cls=aa.Mapper):
+
+                try:
+                    plane_index = (
+                        self.max_log_likelihood_tracer.plane_index_via_redshift_from(
+                            redshift=plane_redshift
+                        )
+                    )
+                    # plane_indexes_with_pixelizations is a list of plane indices that
+                    # have a pixelization, in mapper order. The mapper index for a given
+                    # plane is the position of that plane index within the list, not the
+                    # element at [plane_index].
+                    mapper_index = (
+                        self.max_log_likelihood_tracer.plane_indexes_with_pixelizations.index(
+                            plane_index
+                        )
+                    )
+                except (TypeError, ValueError):
+                    mapper_index = 0
+
+                inversion = self.max_log_likelihood_fit.inversion
+
+                return inversion.max_pixel_centre(mapper_index=mapper_index)
diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py
index 528e2b5cf..edea62c0a 100644
--- a/test_autolens/analysis/test_result.py
+++ b/test_autolens/analysis/test_result.py
@@ -1,385 +1,385 @@
-from pathlib import Path
-import numpy as np
-import pytest
-
-import autofit as af
-import autolens as al
-from autoarray import Array2D
-
-from autolens.analysis import result as res
-from autolens.imaging.model.result import ResultImaging
-
-directory = Path(__file__).resolve().parent
-
-
-def test__max_log_likelihood_tracer(
-    analysis_imaging_7x7,
-    tracer_x2_plane_7x7,
-):
-    model = af.Collection(
-        galaxies=af.Collection(
-            lens=al.Galaxy(redshift=0.5, mass=al.mp.Isothermal),
-            source=al.Galaxy(redshift=1.0, light=al.lp.Sersic),
-        )
-    )
-
-    search = al.m.MockSearch(name="test_search_2")
-
-    result = search.fit(model=model, analysis=analysis_imaging_7x7)
-
-    assert isinstance(result.max_log_likelihood_tracer, al.Tracer)
-    assert isinstance(result.max_log_likelihood_tracer.galaxies[0], al.Galaxy)
-
-
-def test__source_plane_light_profile_centre(analysis_imaging_7x7):
-    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
-
-    source = al.Galaxy(
-        redshift=1.0, light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0)
-    )
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)]
-
-    source_0 = al.Galaxy(
-        redshift=1.0,
-        light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0),
-        light1=al.lp.SersicSph(centre=(3.0, 4.0), intensity=2.0),
-    )
-
-    source_1 = al.Galaxy(
-        redshift=1.0, light=al.lp.SersicSph(centre=(5.0, 6.0), intensity=2.0)
-    )
-
-    tracer = al.Tracer(galaxies=[lens, source_0, source_1])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)]
-
-    source_0 = al.Galaxy(
-        redshift=1.0, light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0)
-    )
-
-    source_1 = al.Galaxy(
-        redshift=2.0, light=al.lp.SersicSph(centre=(5.0, 6.0), intensity=2.0)
-    )
-
-    tracer = al.Tracer(galaxies=[lens, source_0, source_1])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    assert result.source_plane_light_profile_centre_from().in_list == [(5.0, 6.0)]
-
-    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5)])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    assert result.source_plane_light_profile_centre_from() == None
-
-
-def test__source_plane_inversion_centre(analysis_imaging_7x7):
-    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
-
-    pixelization = al.Pixelization(
-        mesh=al.mesh.RectangularUniform((3, 3)),
-        regularization=al.reg.Constant(coefficient=1.0),
-    )
-
-    source = al.Galaxy(redshift=1.0, pixelization=pixelization)
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = ResultImaging(
-        samples_summary=samples_summary, analysis=analysis_imaging_7x7
-    )
-
-    inversion = result.max_log_likelihood_fit.inversion
-
-    assert (
-        result.source_plane_inversion_centre_from().in_list[0]
-        == inversion.max_pixel_centre().in_list[0]
-    )
-
-    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
-    source = al.Galaxy(redshift=1.0)
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = ResultImaging(
-        samples_summary=samples_summary, analysis=analysis_imaging_7x7
-    )
-
-    assert result.source_plane_inversion_centre_from() == None
-
-    lens = al.Galaxy(redshift=0.5, light=al.lp_linear.Sersic())
-    source = al.Galaxy(redshift=1.0)
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = ResultImaging(
-        samples_summary=samples_summary, analysis=analysis_imaging_7x7
-    )
-
-    assert result.source_plane_inversion_centre_from() == None
-
-
-def test__source_plane_centre(analysis_imaging_7x7):
-    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
-
-    pixelization = al.Pixelization(
-        mesh=al.mesh.RectangularUniform((3, 3)),
-        regularization=al.reg.Constant(coefficient=1.0),
-    )
-
-    source = al.Galaxy(
-        redshift=1.0,
-        light=al.lp.SersicSph(centre=(9.0, 8.0), intensity=2.0),
-        pixelization=pixelization,
-    )
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = ResultImaging(
-        samples_summary=samples_summary, analysis=analysis_imaging_7x7
-    )
-
-    assert result.source_plane_centre_from().in_list[0] == pytest.approx(
-        (-0.916666673333, -0.916666), 1.0e-4
-    )
-
-
-def test__image_plane_multiple_image_positions(analysis_imaging_7x7):
-    lens = al.Galaxy(
-        redshift=0.5,
-        mass=al.mp.Isothermal(
-            centre=(0.001, 0.001), einstein_radius=1.0, ell_comps=(0.0, 0.111111)
-        ),
-    )
-
-    source = al.Galaxy(
-        redshift=1.0,
-        light1=al.lp.SersicSph(centre=(0.0, 0.05), intensity=2.0),
-    )
-
-    tracer = al.Tracer(galaxies=[lens, source])
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
-
-    result = ResultImaging(
-        samples_summary=samples_summary, analysis=analysis_imaging_7x7
-    )
-
-    multiple_images = result.image_plane_multiple_image_positions()
-
-    assert pytest.approx((0.968719, 0.366210), 1.0e-2) in multiple_images.in_list
-
-
-def test__positions_threshold_from(analysis_imaging_7x7):
-    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)
-
-    assert result.positions_threshold_from() == pytest.approx(0.0019501455, 1.0e-4)
-    assert result.positions_threshold_from(factor=5.0) == pytest.approx(
-        0.0097507278035, 1.0e-4
-    )
-    assert result.positions_threshold_from(minimum_threshold=10.0) == pytest.approx(
-        10.0, 1.0e-4
-    )
-    assert result.positions_threshold_from(
-        positions=al.Grid2DIrregular([(0.0, 0.0)])
-    ) == pytest.approx(0.0, 1.0e-4)
-
-
-def test__positions_likelihood_from(analysis_imaging_7x7):
-    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)
-
-    positions_likelihood = result.positions_likelihood_from(
-        factor=0.1, minimum_threshold=0.2
-    )
-
-    assert isinstance(positions_likelihood, al.PositionsLH)
-    assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
-
-
-def test__positions_likelihood_from__skip_checks_returns_none_outside_test_mode(
-    monkeypatch, analysis_imaging_7x7,
-):
-    monkeypatch.setenv("PYAUTO_SKIP_CHECKS", "1")
-    monkeypatch.delenv("PYAUTO_TEST_MODE", raising=False)
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=al.Tracer(galaxies=[]))
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    assert result.positions_likelihood_from(factor=0.1, minimum_threshold=0.2) is None
-
-
-def test__positions_likelihood_from__skip_checks_returns_synthetic_in_test_mode(
-    monkeypatch, analysis_imaging_7x7,
-):
-    monkeypatch.setenv("PYAUTO_SKIP_CHECKS", "1")
-    monkeypatch.setenv("PYAUTO_TEST_MODE", "2")
-
-    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=al.Tracer(galaxies=[]))
-    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
-
-    positions_likelihood = result.positions_likelihood_from(
-        factor=0.1, minimum_threshold=0.2
-    )
-
-    assert isinstance(positions_likelihood, al.PositionsLH)
-    assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
-    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__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,
-):
-    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)
-
-    positions_likelihood = result.positions_likelihood_from(
-        factor=0.1, minimum_threshold=0.2, mass_centre_radial_distance_min=0.1
-    )
-
-    assert isinstance(positions_likelihood, al.PositionsLH)
-    assert len(positions_likelihood.positions) == 2
-    assert positions_likelihood.positions[0] == pytest.approx(
-        (-1.00097656e00, 5.63818622e-04), 1.0e-4
-    )
-    assert positions_likelihood.positions[1] == pytest.approx(
-        (1.00097656e00, -5.63818622e-04), 1.0e-4
-    )
-
-
-def test__results_include_mask__available_as_property(
-    analysis_imaging_7x7, masked_imaging_7x7, samples_summary_with_result
-):
-    result = res.ResultDataset(
-        samples_summary=samples_summary_with_result,
-        analysis=analysis_imaging_7x7,
-    )
-
-    assert (result.mask == masked_imaging_7x7.mask).all()
-
-
-def test___image_dict(analysis_imaging_7x7):
-    galaxies = af.ModelInstance()
-    galaxies.lens = al.Galaxy(redshift=0.5)
-    galaxies.source = al.Galaxy(redshift=1.0)
-
-    instance = af.ModelInstance()
-    instance.galaxies = galaxies
-
-    result = ResultImaging(
-        samples_summary=al.m.MockSamplesSummary(max_log_likelihood_instance=instance),
-        analysis=analysis_imaging_7x7,
-    )
-
-    image_dict = result.model_image_galaxy_dict
-
-    assert isinstance(image_dict[str(("galaxies", "lens"))], Array2D)
-    assert isinstance(image_dict[str(("galaxies", "source"))], Array2D)
-
-    result.instance.galaxies.lens = al.Galaxy(redshift=0.5)
-
-    image_dict = result.model_image_galaxy_dict
-
-    assert (image_dict[str(("galaxies", "lens"))].native == np.zeros((7, 7))).all()
-    assert isinstance(image_dict[str(("galaxies", "source"))], Array2D)
+from pathlib import Path
+import numpy as np
+import pytest
+
+import autofit as af
+import autolens as al
+from autoarray import Array2D
+
+from autolens.analysis import result as res
+from autolens.imaging.model.result import ResultImaging
+
+directory = Path(__file__).resolve().parent
+
+
+def test__max_log_likelihood_tracer(
+    analysis_imaging_7x7,
+    tracer_x2_plane_7x7,
+):
+    model = af.Collection(
+        galaxies=af.Collection(
+            lens=al.Galaxy(redshift=0.5, mass=al.mp.Isothermal),
+            source=al.Galaxy(redshift=1.0, light=al.lp.Sersic),
+        )
+    )
+
+    search = al.m.MockSearch(name="test_search_2")
+
+    result = search.fit(model=model, analysis=analysis_imaging_7x7)
+
+    assert isinstance(result.max_log_likelihood_tracer, al.Tracer)
+    assert isinstance(result.max_log_likelihood_tracer.galaxies[0], al.Galaxy)
+
+
+def test__source_plane_light_profile_centre(analysis_imaging_7x7):
+    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
+
+    source = al.Galaxy(
+        redshift=1.0, light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0)
+    )
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)]
+
+    source_0 = al.Galaxy(
+        redshift=1.0,
+        light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0),
+        light1=al.lp.SersicSph(centre=(3.0, 4.0), intensity=2.0),
+    )
+
+    source_1 = al.Galaxy(
+        redshift=1.0, light=al.lp.SersicSph(centre=(5.0, 6.0), intensity=2.0)
+    )
+
+    tracer = al.Tracer(galaxies=[lens, source_0, source_1])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)]
+
+    source_0 = al.Galaxy(
+        redshift=1.0, light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0)
+    )
+
+    source_1 = al.Galaxy(
+        redshift=2.0, light=al.lp.SersicSph(centre=(5.0, 6.0), intensity=2.0)
+    )
+
+    tracer = al.Tracer(galaxies=[lens, source_0, source_1])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    assert result.source_plane_light_profile_centre_from().in_list == [(5.0, 6.0)]
+
+    tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5)])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    assert result.source_plane_light_profile_centre_from() == None
+
+
+def test__source_plane_inversion_centre(analysis_imaging_7x7):
+    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
+
+    pixelization = al.Pixelization(
+        mesh=al.mesh.RectangularUniform((3, 3)),
+        regularization=al.reg.Constant(coefficient=1.0),
+    )
+
+    source = al.Galaxy(redshift=1.0, pixelization=pixelization)
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = ResultImaging(
+        samples_summary=samples_summary, analysis=analysis_imaging_7x7
+    )
+
+    inversion = result.max_log_likelihood_fit.inversion
+
+    assert (
+        result.source_plane_inversion_centre_from().in_list[0]
+        == inversion.max_pixel_centre().in_list[0]
+    )
+
+    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
+    source = al.Galaxy(redshift=1.0)
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = ResultImaging(
+        samples_summary=samples_summary, analysis=analysis_imaging_7x7
+    )
+
+    assert result.source_plane_inversion_centre_from() == None
+
+    lens = al.Galaxy(redshift=0.5, light=al.lp_linear.Sersic())
+    source = al.Galaxy(redshift=1.0)
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = ResultImaging(
+        samples_summary=samples_summary, analysis=analysis_imaging_7x7
+    )
+
+    assert result.source_plane_inversion_centre_from() == None
+
+
+def test__source_plane_centre(analysis_imaging_7x7):
+    lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0))
+
+    pixelization = al.Pixelization(
+        mesh=al.mesh.RectangularUniform((3, 3)),
+        regularization=al.reg.Constant(coefficient=1.0),
+    )
+
+    source = al.Galaxy(
+        redshift=1.0,
+        light=al.lp.SersicSph(centre=(9.0, 8.0), intensity=2.0),
+        pixelization=pixelization,
+    )
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = ResultImaging(
+        samples_summary=samples_summary, analysis=analysis_imaging_7x7
+    )
+
+    assert result.source_plane_centre_from().in_list[0] == pytest.approx(
+        (-0.916666673333, -0.916666), 1.0e-4
+    )
+
+
+def test__image_plane_multiple_image_positions(analysis_imaging_7x7):
+    lens = al.Galaxy(
+        redshift=0.5,
+        mass=al.mp.Isothermal(
+            centre=(0.001, 0.001), einstein_radius=1.0, ell_comps=(0.0, 0.111111)
+        ),
+    )
+
+    source = al.Galaxy(
+        redshift=1.0,
+        light1=al.lp.SersicSph(centre=(0.0, 0.05), intensity=2.0),
+    )
+
+    tracer = al.Tracer(galaxies=[lens, source])
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = ResultImaging(
+        samples_summary=samples_summary, analysis=analysis_imaging_7x7
+    )
+
+    multiple_images = result.image_plane_multiple_image_positions()
+
+    assert pytest.approx((0.968719, 0.366210), 1.0e-2) in multiple_images.in_list
+
+
+def test__positions_threshold_from(analysis_imaging_7x7):
+    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)
+
+    assert result.positions_threshold_from() == pytest.approx(0.0019501455, 1.0e-4)
+    assert result.positions_threshold_from(factor=5.0) == pytest.approx(
+        0.0097507278035, 1.0e-4
+    )
+    assert result.positions_threshold_from(minimum_threshold=10.0) == pytest.approx(
+        10.0, 1.0e-4
+    )
+    assert result.positions_threshold_from(
+        positions=al.Grid2DIrregular([(0.0, 0.0)])
+    ) == pytest.approx(0.0, 1.0e-4)
+
+
+def test__positions_likelihood_from(analysis_imaging_7x7):
+    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)
+
+    positions_likelihood = result.positions_likelihood_from(
+        factor=0.1, minimum_threshold=0.2
+    )
+
+    assert isinstance(positions_likelihood, al.PositionsLH)
+    assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
+
+
+def test__positions_likelihood_from__skip_checks_returns_none_outside_test_mode(
+    monkeypatch, analysis_imaging_7x7,
+):
+    monkeypatch.setenv("PYAUTO_SKIP_CHECKS", "1")
+    monkeypatch.delenv("PYAUTO_TEST_MODE", raising=False)
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=al.Tracer(galaxies=[]))
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    assert result.positions_likelihood_from(factor=0.1, minimum_threshold=0.2) is None
+
+
+def test__positions_likelihood_from__skip_checks_returns_synthetic_in_test_mode(
+    monkeypatch, analysis_imaging_7x7,
+):
+    monkeypatch.setenv("PYAUTO_SKIP_CHECKS", "1")
+    monkeypatch.setenv("PYAUTO_TEST_MODE", "2")
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=al.Tracer(galaxies=[]))
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    positions_likelihood = result.positions_likelihood_from(
+        factor=0.1, minimum_threshold=0.2
+    )
+
+    assert isinstance(positions_likelihood, al.PositionsLH)
+    assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
+    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__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,
+):
+    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)
+
+    positions_likelihood = result.positions_likelihood_from(
+        factor=0.1, minimum_threshold=0.2, mass_centre_radial_distance_min=0.1
+    )
+
+    assert isinstance(positions_likelihood, al.PositionsLH)
+    assert len(positions_likelihood.positions) == 2
+    assert positions_likelihood.positions[0] == pytest.approx(
+        (-1.00097656e00, 5.63818622e-04), 1.0e-4
+    )
+    assert positions_likelihood.positions[1] == pytest.approx(
+        (1.00097656e00, -5.63818622e-04), 1.0e-4
+    )
+
+
+def test__results_include_mask__available_as_property(
+    analysis_imaging_7x7, masked_imaging_7x7, samples_summary_with_result
+):
+    result = res.ResultDataset(
+        samples_summary=samples_summary_with_result,
+        analysis=analysis_imaging_7x7,
+    )
+
+    assert (result.mask == masked_imaging_7x7.mask).all()
+
+
+def test___image_dict(analysis_imaging_7x7):
+    galaxies = af.ModelInstance()
+    galaxies.lens = al.Galaxy(redshift=0.5)
+    galaxies.source = al.Galaxy(redshift=1.0)
+
+    instance = af.ModelInstance()
+    instance.galaxies = galaxies
+
+    result = ResultImaging(
+        samples_summary=al.m.MockSamplesSummary(max_log_likelihood_instance=instance),
+        analysis=analysis_imaging_7x7,
+    )
+
+    image_dict = result.model_image_galaxy_dict
+
+    assert isinstance(image_dict[str(("galaxies", "lens"))], Array2D)
+    assert isinstance(image_dict[str(("galaxies", "source"))], Array2D)
+
+    result.instance.galaxies.lens = al.Galaxy(redshift=0.5)
+
+    image_dict = result.model_image_galaxy_dict
+
+    assert (image_dict[str(("galaxies", "lens"))].native == np.zeros((7, 7))).all()
+    assert isinstance(image_dict[str(("galaxies", "source"))], Array2D)
 
 
 def test__positions_likelihood_from__loads_cached_positions_on_second_call(
@@ -389,6 +389,9 @@ class _StubPaths:
         def __init__(self, files_path):
             self._files_path = files_path
 
+        def preserve_in_zip(self, file_path):
+            """No zip in the stub — mirrors AbstractPaths' no-zip no-op."""
+
     tracer = al.Tracer(
         galaxies=[
             al.Galaxy(