diff --git a/autolens/__init__.py b/autolens/__init__.py
index c1a45fd86..94a49c0a3 100644
--- a/autolens/__init__.py
+++ b/autolens/__init__.py
@@ -123,4 +123,4 @@
 
 conf.instance.register(__file__)
 
-__version__ = "2024.11.6.1"
+__version__ = "2024.11.13.2"
diff --git a/autolens/aggregator/tracer.py b/autolens/aggregator/tracer.py
index e465c00ec..a0d1efb3a 100644
--- a/autolens/aggregator/tracer.py
+++ b/autolens/aggregator/tracer.py
@@ -58,17 +58,18 @@ def _tracer_from(
 
     tracer = Tracer(galaxies=galaxies, cosmology=cosmology)
 
-    if len(fit.children) > 0:
-        logger.info(
-            """
-            Using database for a fit with multiple summed Analysis objects.
-
-            Tracer objects do not fully support this yet (e.g. model parameters which vary over analyses may be incorrect)
-            so proceed with caution!
-            """
-        )
-
-        return [tracer] * len(fit.children)
+    if fit.children is not None:
+        if len(fit.children) > 0:
+            logger.info(
+                """
+                Using database for a fit with multiple summed Analysis objects.
+    
+                Tracer objects do not fully support this yet (e.g. model parameters which vary over analyses may be incorrect)
+                so proceed with caution!
+                """
+            )
+
+            return [tracer] * len(fit.children)
 
     return [tracer]
 
diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py
index 3ec792ec8..fa9ed4d11 100644
--- a/autolens/analysis/result.py
+++ b/autolens/analysis/result.py
@@ -1,3 +1,4 @@
+import logging
 import os
 import numpy as np
 from typing import Optional, Union
@@ -15,6 +16,8 @@
 from autolens.lens.tracer import Tracer
 from autolens.point.solver import PointSolver
 
+logger = logging.getLogger(__name__)
+
 
 class Result(AgResultDataset):
     @property
@@ -101,8 +104,76 @@ def image_plane_multiple_image_positions(self) -> aa.Grid2DIrregular:
             source_plane_coordinate=self.source_plane_centre.in_list[0],
         )
 
+        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, 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
+        ----------
+        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.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),
+            )
+
+            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,
@@ -168,6 +239,7 @@ def positions_likelihood_from(
         minimum_threshold=None,
         use_resample=False,
         positions: Optional[aa.Grid2DIrregular] = None,
+        mass_centre_radial_distance_min: float = None,
     ) -> Union[PositionsLHPenalty, PositionsLHResample]:
         """
         Returns a `PositionsLH` object from the result of a lens model-fit, where the maximum log likelihood mass
@@ -178,6 +250,11 @@ def positions_likelihood_from(
         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
@@ -193,6 +270,10 @@ def positions_likelihood_from(
         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.
 
         Returns
         -------
@@ -207,6 +288,18 @@ def positions_likelihood_from(
             if positions is None
             else positions
         )
+
+        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]
+
         threshold = self.positions_threshold_from(
             factor=factor, minimum_threshold=minimum_threshold, positions=positions
         )
diff --git a/autolens/config/visualize/plots.yaml b/autolens/config/visualize/plots.yaml
index d1c85c0ea..372d27db6 100644
--- a/autolens/config/visualize/plots.yaml
+++ b/autolens/config/visualize/plots.yaml
@@ -55,7 +55,7 @@
     all_at_end_png: true                     # Plot all individual plots listed below as .png (even if False)?
     all_at_end_fits: true                    # Plot all individual plots listed below as .fits (even if False)?
     all_at_end_pdf: false                    # Plot all individual plots listed below as publication-quality .pdf (even if False)?
-    errors: false
+    reconstruction_noise_map: false
     reconstructed_image: false
     reconstruction: false
     regularization_weights: false
diff --git a/autolens/imaging/plot/fit_imaging_plotters.py b/autolens/imaging/plot/fit_imaging_plotters.py
index 463adfaa2..5ddfb2295 100644
--- a/autolens/imaging/plot/fit_imaging_plotters.py
+++ b/autolens/imaging/plot/fit_imaging_plotters.py
@@ -148,7 +148,7 @@ def figures_2d_of_planes(
         subtracted_image: bool = False,
         model_image: bool = False,
         plane_image: bool = False,
-        plane_errors: bool = False,
+        plane_noise_map: bool = False,
         plane_signal_to_noise_map: bool = False,
         use_source_vmax: bool = False,
         zoom_to_brightest: bool = True,
@@ -180,12 +180,12 @@ def figures_2d_of_planes(
             Whether to make a 2D plot (via `imshow`) of the image of a plane in its source-plane (e.g. unlensed).
             Depending on how the fit is performed, this could either be an image of light profiles of the reconstruction
             of an `Inversion`.
-        plane_errors
-            Whether to make a 2D plot (via `imshow`) of the errors of a plane in its source-plane, where the
-            errors can only be computed when a pixelized source reconstruction is performed and they correspond to
-            the errors in each reconstructed pixel as given by the inverse curvature matrix.
+        plane_noise_map
+            Whether to make a 2D plot of the noise-map of a plane in its source-plane, where the
+            noise map can only be computed when a pixelized source reconstruction is performed and they correspond to
+            the noise map in each reconstructed pixel as given by the inverse curvature matrix.
         plane_signal_to_noise_map
-            Whether to make a 2D plot (via `imshow`) of the signal-to-noise map of a plane in its source-plane,
+            Whether to make a 2D plot of the signal-to-noise map of a plane in its source-plane,
             where the signal-to-noise map values can only be computed when a pixelized source reconstruction and they
             are the ratio of reconstructed flux to error in each pixel.
         use_source_vmax
@@ -295,12 +295,6 @@ def figures_2d_of_planes(
 
                     pix = self.tracer.planes[plane_index].cls_list_from(cls=aa.Pixelization)[0]
 
-                    if isinstance(pix.mesh, aa.mesh.Delaunay):
-                        try:
-                            self.mat_plot_2d.cmap.kwargs.pop("vmax")
-                        except KeyError:
-                            pass
-
                     inversion_plotter = self.inversion_plotter_of_plane(
                         plane_index=plane_index
                     )
@@ -318,7 +312,7 @@ def figures_2d_of_planes(
                 except KeyError:
                     pass
 
-            if plane_errors:
+            if plane_noise_map:
 
                 if self.tracer.planes[plane_index].has(cls=aa.Pixelization):
 
@@ -328,7 +322,7 @@ def figures_2d_of_planes(
 
                     inversion_plotter.figures_2d_of_pixelization(
                         pixelization_index=0,
-                        errors=True,
+                        reconstruction_noise_map=True,
                         zoom_to_brightest=zoom_to_brightest,
                         interpolate_to_uniform=interpolate_to_uniform
                     )
diff --git a/autolens/imaging/simulator.py b/autolens/imaging/simulator.py
index 2bc49bfeb..e740d254a 100644
--- a/autolens/imaging/simulator.py
+++ b/autolens/imaging/simulator.py
@@ -9,21 +9,30 @@ class SimulatorImaging(aa.SimulatorImaging):
 
     def via_tracer_from(self, tracer : Tracer, grid : aa.type.Grid2DLike) -> aa.Imaging:
         """
-        Simulate an `Imaging` dataset from an input tracer and grid.
+        Simulate an `Imaging` dataset from an input `Tracer` object and a 2D grid of (y,x) coordinates.
 
-        The tracer is used to perform ray-tracing and generate the image of the strong lens galaxies (e.g.
-        the lens light, lensed source light, etc) which is simulated.
+        The mass profiles of each galaxy in the tracer are used to perform ray-tracing of the input 2D grid and
+        their light profiles are used to generate the image of the galaxies which are simulated.
 
         The steps of the `SimulatorImaging` simulation process (e.g. PSF convolution, noise addition) are
-        described in the `SimulatorImaging` `__init__` method docstring.
+        described in the `SimulatorImaging` `__init__` method docstring, found in the PyAutoArray project.
+
+        If one of more galaxy light profiles are a `LightProfileSNR` object, the `intensity` of the light profile is
+        automatically set such that the signal-to-noise ratio of the light profile is equal to its input
+        `signal_to_noise_ratio` value.
+
+        For example, if a `LightProfileSNR` object has a `signal_to_noise_ratio` of 5.0, the intensity of the light
+        profile is set such that the peak surface brightness of the profile is 5.0 times the background noise level of
+        the image.
 
         Parameters
         ----------
         tracer
             The tracer, which describes the ray-tracing and strong lens configuration used to simulate the imaging
-            dataset.
+            dataset as well as the light profiles of the galaxies used to simulate the image of the galaxies.
         grid
-            The image-plane grid which the image of the strong lens is generated on.
+            The 2D grid of (y,x) coordinates which the mass profiles of the galaxies in the tracer are ray-traced using
+            in order to generate the image of the galaxies via their light profiles.
         """
 
         tracer.set_snr_of_snr_light_profiles(
@@ -44,21 +53,30 @@ def via_tracer_from(self, tracer : Tracer, grid : aa.type.Grid2DLike) -> aa.Imag
 
     def via_galaxies_from(self, galaxies : List[ag.Galaxy], grid : aa.type.Grid2DLike) -> aa.Imaging:
         """
-        Simulate an `Imaging` dataset from an input list of galaxies and grid.
+        Simulate an `Imaging` dataset from an input list of `Galaxy` objects and a 2D grid of (y,x) coordinates.
 
-        The galaxies are used to create a tracer, which performs ray-tracing and generate the image of the strong lens
-        galaxies (e.g. the lens light, lensed source light, etc) which is simulated.
+        The galaxies are used to create a `Tracer`. The mass profiles of each galaxy in the tracer are used to
+        perform ray-tracing of the input 2D grid and their light profiles are used to generate the image of the
+        galaxies which are simulated.
 
         The steps of the `SimulatorImaging` simulation process (e.g. PSF convolution, noise addition) are
         described in the `SimulatorImaging` `__init__` method docstring.
 
+        If one of more galaxy light profiles are a `LightProfileSNR` object, the `intensity` of the light profile is
+        automatically set such that the signal-to-noise ratio of the light profile is equal to its input
+        `signal_to_noise_ratio` value.
+
+        For example, if a `LightProfileSNR` object has a `signal_to_noise_ratio` of 5.0, the intensity of the light
+        profile is set such that the peak surface brightness of the profile is 5.0 times the background noise level of
+        the image.
+
         Parameters
         ----------
         galaxies
             The galaxies used to create the tracer, which describes the ray-tracing and strong lens configuration
             used to simulate the imaging dataset.
         grid
-            The image-plane grid which the image of the strong lens is generated on.
+            The image-plane 2D grid of (y,x) coordinates grid which the image of the strong lens is generated on.
         """
 
         tracer = Tracer(galaxies=galaxies)
@@ -87,7 +105,7 @@ def via_deflections_and_galaxies_from(self, deflections : aa.VectorYX2D, galaxie
             The galaxies used to create the tracer, which describes the ray-tracing and strong lens configuration
             used to simulate the imaging dataset.
         grid
-            The image-plane grid which the image of the strong lens is generated on.
+            The image-plane 2D grid of (y,x) coordinates grid which the image of the strong lens is generated on.
         """
         grid = aa.Grid2D.uniform(
             shape_native=deflections.shape_native,
@@ -105,10 +123,10 @@ def via_source_image_from(self, tracer : Tracer, grid : aa.type.Grid2DLike, sour
         Simulate an `Imaging` dataset from an input image of a source galaxy.
 
         This input image is on a uniform and regular 2D array, meaning it can simulate the source's irregular
-        and assymetric source galaxy morphological features.
+        and asymmetric source galaxy morphological features.
 
         The typical use case is inputting the image of an irregular galaxy in the source-plane (whose values are
-        on a uniform array) and using this function computing the lensed image of this source galaxy.
+        on a uniform array) and using this function to compute the lensed image of this source galaxy.
 
         The tracer is used to perform ray-tracing and generate the image of the strong lens galaxies (e.g.
         the lens light, lensed source light, etc) which is simulated.
@@ -125,7 +143,7 @@ def via_source_image_from(self, tracer : Tracer, grid : aa.type.Grid2DLike, sour
             The tracer, which describes the ray-tracing and strong lens configuration used to simulate the imaging
             dataset.
         grid
-            The image-plane grid which the image of the strong lens is generated on.
+            The image-plane 2D grid of (y,x) coordinates grid which the image of the strong lens is generated on.
         source_image
             The image of the source-plane and source galaxy which is interpolated to compute the lensed image.
         """
diff --git a/autolens/interferometer/plot/fit_interferometer_plotters.py b/autolens/interferometer/plot/fit_interferometer_plotters.py
index 0bfd1fb7e..73cc95986 100644
--- a/autolens/interferometer/plot/fit_interferometer_plotters.py
+++ b/autolens/interferometer/plot/fit_interferometer_plotters.py
@@ -375,7 +375,7 @@ def figures_2d_of_planes(
         self,
         plane_index: Optional[int] = None,
         plane_image: bool = False,
-        plane_errors: bool = False,
+        plane_noise_map: bool = False,
         plane_signal_to_noise_map: bool = False,
         zoom_to_brightest: bool = True,
         interpolate_to_uniform: bool = False,
@@ -398,12 +398,12 @@ def figures_2d_of_planes(
             Whether to make a 2D plot (via `imshow`) of the image of a plane in its source-plane (e.g. unlensed).
             Depending on how the fit is performed, this could either be an image of light profiles of the reconstruction
             of an `Inversion`.
-        plane_errors
-            Whether to make a 2D plot (via `imshow`) of the errors of a plane in its source-plane, where the
-            errors can only be computed when a pixelized source reconstruction is performed and they correspond to
-            the errors in each reconstructed pixel as given by the inverse curvature matrix.
+        plane_noise_map
+            Whether to make a 2D plot of the noise-map of a plane in its source-plane, where the
+            noise map can only be computed when a pixelized source reconstruction is performed and they correspond to
+            the noise map in each reconstructed pixel as given by the inverse curvature matrix.
         plane_signal_to_noise_map
-            Whether to make a 2D plot (via `imshow`) of the signal-to-noise map of a plane in its source-plane,
+            Whether to make a 2D plot of the signal-to-noise map of a plane in its source-plane,
             where the signal-to-noise map values can only be computed when a pixelized source reconstruction and they
             are the ratio of reconstructed flux to error in each pixel.
         zoom_to_brightest
@@ -430,7 +430,7 @@ def figures_2d_of_planes(
                     interpolate_to_uniform=interpolate_to_uniform,
                 )
 
-        if plane_errors:
+        if plane_noise_map:
             if self.tracer.planes[plane_index].has(cls=aa.Pixelization):
                 inversion_plotter = self.inversion_plotter_of_plane(
                     plane_index=plane_index
@@ -438,7 +438,7 @@ def figures_2d_of_planes(
 
                 inversion_plotter.figures_2d_of_pixelization(
                     pixelization_index=0,
-                    errors=True,
+                    reconstruction_noise_map=True,
                     zoom_to_brightest=zoom_to_brightest,
                     interpolate_to_uniform=interpolate_to_uniform,
                 )
diff --git a/autolens/interferometer/simulator.py b/autolens/interferometer/simulator.py
index b7592abbb..1245bb134 100644
--- a/autolens/interferometer/simulator.py
+++ b/autolens/interferometer/simulator.py
@@ -18,7 +18,7 @@ def via_tracer_from(self, tracer, grid):
             An arrays representing the effective exposure time of each pixel.
         psf: PSF
             An arrays describing the PSF the simulated image is blurred with.
-        add_poisson_noise: Bool
+        add_poisson_noise_to_data: Bool
             If `True` poisson noise_maps is simulated and added to the image, based on the total counts in each image
             pixel
         noise_seed: int
diff --git a/autolens/lens/sensitivity.py b/autolens/lens/sensitivity.py
index 7c969aa20..6511ba858 100644
--- a/autolens/lens/sensitivity.py
+++ b/autolens/lens/sensitivity.py
@@ -361,7 +361,6 @@ def set_auto_filename(
         return False
 
     def subplot_sensitivity(self):
-
         log_likelihoods = self.result.figure_of_merit_array(
             use_log_evidences=False,
             remove_zeros=True,
@@ -375,7 +374,7 @@ def subplot_sensitivity(self):
         except TypeError:
             log_evidences = np.zeros_like(log_likelihoods)
 
-        self.open_subplot_figure(number_subplots=8, subplot_shape=(2,4))
+        self.open_subplot_figure(number_subplots=8, subplot_shape=(2, 4))
 
         plotter = aplt.Array2DPlotter(
             array=self.data_subtracted,
@@ -398,10 +397,7 @@ def subplot_sensitivity(self):
 
         above_threshold = np.where(log_likelihoods > 5.0, 1.0, 0.0)
 
-        above_threshold = aa.Array2D(
-            values=above_threshold,
-            mask=log_likelihoods.mask
-        )
+        above_threshold = aa.Array2D(values=above_threshold, mask=log_likelihoods.mask)
 
         self.mat_plot_2d.plot_array(
             array=above_threshold,
@@ -410,16 +406,32 @@ def subplot_sensitivity(self):
         )
 
         try:
-            log_evidences_base = self.result._array_2d_from(self.result.log_evidences_base)
-            log_evidences_perturbed = self.result._array_2d_from(self.result.log_evidences_perturbed)
+            log_evidences_base = self.result._array_2d_from(
+                self.result.log_evidences_base
+            )
+            log_evidences_perturbed = self.result._array_2d_from(
+                self.result.log_evidences_perturbed
+            )
 
-            log_evidences_base_min = np.nanmin(np.where(log_evidences_base == 0, np.nan, log_evidences_base))
-            log_evidences_base_max = np.nanmax(np.where(log_evidences_base == 0, np.nan, log_evidences_base))
-            log_evidences_perturbed_min = np.nanmin(np.where(log_evidences_perturbed == 0, np.nan, log_evidences_perturbed))
-            log_evidences_perturbed_max = np.nanmax(np.where(log_evidences_perturbed == 0, np.nan, log_evidences_perturbed))
+            log_evidences_base_min = np.nanmin(
+                np.where(log_evidences_base == 0, np.nan, log_evidences_base)
+            )
+            log_evidences_base_max = np.nanmax(
+                np.where(log_evidences_base == 0, np.nan, log_evidences_base)
+            )
+            log_evidences_perturbed_min = np.nanmin(
+                np.where(log_evidences_perturbed == 0, np.nan, log_evidences_perturbed)
+            )
+            log_evidences_perturbed_max = np.nanmax(
+                np.where(log_evidences_perturbed == 0, np.nan, log_evidences_perturbed)
+            )
 
-            self.mat_plot_2d.cmap.kwargs["vmin"] = np.min([log_evidences_base_min, log_evidences_perturbed_min])
-            self.mat_plot_2d.cmap.kwargs["vmax"] = np.max([log_evidences_base_max, log_evidences_perturbed_max])
+            self.mat_plot_2d.cmap.kwargs["vmin"] = np.min(
+                [log_evidences_base_min, log_evidences_perturbed_min]
+            )
+            self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(
+                [log_evidences_base_max, log_evidences_perturbed_max]
+            )
 
             self.mat_plot_2d.plot_array(
                 array=log_evidences_base,
@@ -431,21 +443,36 @@ def subplot_sensitivity(self):
                 array=log_evidences_perturbed,
                 visuals_2d=self.visuals_2d,
                 auto_labels=AutoLabels(title="Log Evidence Perturb"),
-
             )
         except TypeError:
             pass
-        
-        log_likelihoods_base = self.result._array_2d_from(self.result.log_likelihoods_base)
-        log_likelihoods_perturbed = self.result._array_2d_from(self.result.log_likelihoods_perturbed)
 
-        log_likelihoods_base_min = np.nanmin(np.where(log_likelihoods_base == 0, np.nan, log_likelihoods_base))
-        log_likelihoods_base_max = np.nanmax(np.where(log_likelihoods_base == 0, np.nan, log_likelihoods_base))
-        log_likelihoods_perturbed_min = np.nanmin(np.where(log_likelihoods_perturbed == 0, np.nan, log_likelihoods_perturbed))
-        log_likelihoods_perturbed_max = np.nanmax(np.where(log_likelihoods_perturbed == 0, np.nan, log_likelihoods_perturbed))
+        log_likelihoods_base = self.result._array_2d_from(
+            self.result.log_likelihoods_base
+        )
+        log_likelihoods_perturbed = self.result._array_2d_from(
+            self.result.log_likelihoods_perturbed
+        )
+
+        log_likelihoods_base_min = np.nanmin(
+            np.where(log_likelihoods_base == 0, np.nan, log_likelihoods_base)
+        )
+        log_likelihoods_base_max = np.nanmax(
+            np.where(log_likelihoods_base == 0, np.nan, log_likelihoods_base)
+        )
+        log_likelihoods_perturbed_min = np.nanmin(
+            np.where(log_likelihoods_perturbed == 0, np.nan, log_likelihoods_perturbed)
+        )
+        log_likelihoods_perturbed_max = np.nanmax(
+            np.where(log_likelihoods_perturbed == 0, np.nan, log_likelihoods_perturbed)
+        )
 
-        self.mat_plot_2d.cmap.kwargs["vmin"] = np.min([log_likelihoods_base_min, log_likelihoods_perturbed_min])
-        self.mat_plot_2d.cmap.kwargs["vmax"] = np.max([log_likelihoods_base_max, log_likelihoods_perturbed_max])
+        self.mat_plot_2d.cmap.kwargs["vmin"] = np.min(
+            [log_likelihoods_base_min, log_likelihoods_perturbed_min]
+        )
+        self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(
+            [log_likelihoods_base_max, log_likelihoods_perturbed_max]
+        )
 
         self.mat_plot_2d.plot_array(
             array=log_likelihoods_base,
diff --git a/autolens/plot/__init__.py b/autolens/plot/__init__.py
index 94b640250..57dd8f2a1 100644
--- a/autolens/plot/__init__.py
+++ b/autolens/plot/__init__.py
@@ -28,6 +28,7 @@
     GridPlot,
     VectorYXQuiver,
     PatchOverlay,
+    DelaunayDrawer,
     VoronoiDrawer,
     OriginScatter,
     MaskScatter,
diff --git a/autolens/point/dataset.py b/autolens/point/dataset.py
index 3214e70f7..75e084fb1 100644
--- a/autolens/point/dataset.py
+++ b/autolens/point/dataset.py
@@ -84,11 +84,10 @@ def info(self) -> str:
         info += f"fluxes_noise_map : {self.fluxes_noise_map}\n"
         return info
 
-    def extent_from(self, buffer : float = 0.1):
-
+    def extent_from(self, buffer: float = 0.1):
         y_max = max(self.positions[:, 0]) + buffer
         y_min = min(self.positions[:, 0]) - buffer
         x_max = max(self.positions[:, 1]) + buffer
         x_min = min(self.positions[:, 1]) - buffer
 
-        return [y_min, y_max, x_min, x_max]
\ No newline at end of file
+        return [y_min, y_max, x_min, x_max]
diff --git a/autolens/point/model/analysis.py b/autolens/point/model/analysis.py
index 3dd23e6e1..b5731f0d0 100644
--- a/autolens/point/model/analysis.py
+++ b/autolens/point/model/analysis.py
@@ -24,7 +24,6 @@
 
 
 class AnalysisPoint(AgAnalysis, AnalysisLens):
-
     Visualizer = VisualizerPoint
     Result = ResultPoint
 
diff --git a/autolens/point/model/plotter_interface.py b/autolens/point/model/plotter_interface.py
index a1962fd83..51417f66a 100644
--- a/autolens/point/model/plotter_interface.py
+++ b/autolens/point/model/plotter_interface.py
@@ -11,7 +11,6 @@
 
 
 class PlotterInterfacePoint(PlotterInterface):
-
     def dataset_point(self, dataset: PointDataset):
         """
         Output visualization of an `PointDataset` dataset, typically before a model-fit is performed.
@@ -54,7 +53,10 @@ def should_plot(name):
             dataset_plotter.subplot_dataset()
 
     def fit_point(
-        self, fit: FitPointDataset, during_analysis: bool, subfolders: str = "fit_dataset"
+        self,
+        fit: FitPointDataset,
+        during_analysis: bool,
+        subfolders: str = "fit_dataset",
     ):
         """
         Visualizes a `FitPointDataset` object, which fits an imaging dataset.
@@ -104,7 +106,6 @@ def should_plot(name):
             fit_plotter.subplot_fit()
 
         if not during_analysis and should_plot("all_at_end_png"):
-
             mat_plot_2d = self.mat_plot_2d_from(
                 subfolders=path.join("fit_dataset", "end"),
             )
@@ -113,7 +114,4 @@ def should_plot(name):
                 fit=fit, mat_plot_2d=mat_plot_2d, include_2d=self.include_2d
             )
 
-            fit_plotter.figures_2d(
-                positions=True,
-                fluxes=True
-            )
+            fit_plotter.figures_2d(positions=True, fluxes=True)
diff --git a/autolens/point/model/visualizer.py b/autolens/point/model/visualizer.py
index b19ea343c..fe86e98f8 100644
--- a/autolens/point/model/visualizer.py
+++ b/autolens/point/model/visualizer.py
@@ -5,7 +5,6 @@
 
 
 class VisualizerPoint(af.Visualizer):
-
     @staticmethod
     def visualize_before_fit(
         analysis,
@@ -76,7 +75,9 @@ def visualize(
 
         tracer = fit.tracer
 
-        grid = ag.Grid2D.from_extent(extent=fit.dataset.extent_from(), shape_native=(100, 100))
+        grid = ag.Grid2D.from_extent(
+            extent=fit.dataset.extent_from(), shape_native=(100, 100)
+        )
 
         plotter_interface.tracer(
             tracer=tracer, grid=grid, during_analysis=during_analysis
diff --git a/autolens/point/plot/fit_point_plotters.py b/autolens/point/plot/fit_point_plotters.py
index 3ac8bf299..8876bcd84 100644
--- a/autolens/point/plot/fit_point_plotters.py
+++ b/autolens/point/plot/fit_point_plotters.py
@@ -41,25 +41,36 @@ def figures_2d(self, positions: bool = False, fluxes: bool = False):
             )
 
             if self.mat_plot_2d.axis.kwargs.get("extent") is None:
-
                 buffer = 0.1
 
-                y_max = max(
-                    max(self.fit.dataset.positions[:, 0]),
-                    max(self.fit.positions.model_data[:, 0]),
-                ) + buffer
-                y_min = min(
-                    min(self.fit.dataset.positions[:, 0]),
-                    min(self.fit.positions.model_data[:, 0]),
-                ) - buffer
-                x_max = max(
-                    max(self.fit.dataset.positions[:, 1]),
-                    max(self.fit.positions.model_data[:, 1]),
-                ) + buffer
-                x_min = min(
-                    min(self.fit.dataset.positions[:, 1]),
-                    min(self.fit.positions.model_data[:, 1]),
-                ) - buffer
+                y_max = (
+                    max(
+                        max(self.fit.dataset.positions[:, 0]),
+                        max(self.fit.positions.model_data[:, 0]),
+                    )
+                    + buffer
+                )
+                y_min = (
+                    min(
+                        min(self.fit.dataset.positions[:, 0]),
+                        min(self.fit.positions.model_data[:, 0]),
+                    )
+                    - buffer
+                )
+                x_max = (
+                    max(
+                        max(self.fit.dataset.positions[:, 1]),
+                        max(self.fit.positions.model_data[:, 1]),
+                    )
+                    + buffer
+                )
+                x_min = (
+                    min(
+                        min(self.fit.dataset.positions[:, 1]),
+                        min(self.fit.positions.model_data[:, 1]),
+                    )
+                    - buffer
+                )
 
                 extent = [y_min, y_max, x_min, x_max]
 
diff --git a/autolens/point/solver/shape_solver.py b/autolens/point/solver/shape_solver.py
index 7c8a040e9..f85043add 100644
--- a/autolens/point/solver/shape_solver.py
+++ b/autolens/point/solver/shape_solver.py
@@ -4,24 +4,27 @@
 from typing import Tuple, List, Iterator, Type, Optional
 
 import autoarray as aa
-from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles
 
 from autoarray.structures.triangles.shape import Shape
 from autofit.jax_wrapper import jit, use_jax, numpy as np, register_pytree_node_class
 
 try:
     if use_jax:
-        from autoarray.structures.triangles.jax_array import (
-            ArrayTriangles,
+        from autoarray.structures.triangles.coordinate_array.jax_coordinate_array import (
+            CoordinateArrayTriangles,
             MAX_CONTAINING_SIZE,
         )
     else:
-        from autoarray.structures.triangles.array import ArrayTriangles
+        from autoarray.structures.triangles.coordinate_array.coordinate_array import (
+            CoordinateArrayTriangles,
+        )
 
         MAX_CONTAINING_SIZE = None
 
 except ImportError:
-    from autoarray.structures.triangles.array import ArrayTriangles
+    from autoarray.structures.triangles.coordinate_array.coordinate_array import (
+        CoordinateArrayTriangles,
+    )
 
     MAX_CONTAINING_SIZE = None
 
@@ -74,6 +77,7 @@ def for_grid(
         magnification_threshold=0.1,
         array_triangles_cls: Type[AbstractTriangles] = CoordinateArrayTriangles,
         max_containing_size=MAX_CONTAINING_SIZE,
+        neighbor_degree: int = 1,
     ):
         """
         Create a solver for a given grid.
@@ -94,6 +98,8 @@ def for_grid(
         max_containing_size
             Only applies to JAX. This is the maximum number of multiple images expected.
             We need to know this in advance to allocate memory for the JAX array.
+        neighbor_degree
+            The number of times recursively add neighbors for the triangles that contain
 
         Returns
         -------
@@ -109,6 +115,64 @@ def for_grid(
         x_min = x.min()
         x_max = x.max()
 
+        return cls.for_limits_and_scale(
+            y_min=y_min,
+            y_max=y_max,
+            x_min=x_min,
+            x_max=x_max,
+            scale=scale,
+            pixel_scale_precision=pixel_scale_precision,
+            magnification_threshold=magnification_threshold,
+            array_triangles_cls=array_triangles_cls,
+            max_containing_size=max_containing_size,
+            neighbor_degree=neighbor_degree,
+        )
+
+    @classmethod
+    def for_limits_and_scale(
+        cls,
+        y_min=-1.0,
+        y_max=1.0,
+        x_min=-1.0,
+        x_max=1.0,
+        scale=0.1,
+        pixel_scale_precision: float = 0.001,
+        magnification_threshold=0.1,
+        array_triangles_cls: Type[AbstractTriangles] = CoordinateArrayTriangles,
+        max_containing_size=MAX_CONTAINING_SIZE,
+        neighbor_degree: int = 1,
+    ):
+        """
+        Create a solver for a given grid.
+
+        The grid defines the limits of the image plane and the pixel scale.
+
+        Parameters
+        ----------
+        y_min
+        y_max
+        x_min
+        x_max
+            The limits of the image plane in pixels.
+        scale
+            The pixel scale of the image plane. The initial triangles have this side length.
+        pixel_scale_precision
+            The precision to which the triangles should be subdivided.
+        magnification_threshold
+            The threshold for the magnification under which multiple images are filtered.
+        array_triangles_cls
+            The class to use for the triangles. JAX is used implicitly if USE_JAX=1 and
+            jax is installed.
+        max_containing_size
+            Only applies to JAX. This is the maximum number of multiple images expected.
+            We need to know this in advance to allocate memory for the JAX array.
+        neighbor_degree
+            The number of times recursively add neighbors for the triangles that contain
+
+        Returns
+        -------
+        The solver.
+        """
         initial_triangles = array_triangles_cls.for_limits_and_scale(
             y_min=y_min,
             y_max=y_max,
@@ -123,6 +187,7 @@ def for_grid(
             initial_triangles=initial_triangles,
             pixel_scale_precision=pixel_scale_precision,
             magnification_threshold=magnification_threshold,
+            neighbor_degree=neighbor_degree,
         )
 
     @property
@@ -233,12 +298,11 @@ def _filter_low_magnification(
         mask = np.abs(magnifications.array) > self.magnification_threshold
         return np.where(mask[:, None], points, np.nan)
 
-    def _filtered_triangles(
+    def _source_triangles(
         self,
         tracer: OperateDeflections,
         triangles: aa.AbstractTriangles,
         source_plane_redshift,
-        shape: Shape,
     ):
         """
         Filter the triangles to keep only those that meet the solver condition
@@ -248,11 +312,7 @@ def _filtered_triangles(
             grid=aa.Grid2DIrregular(triangles.vertices),
             source_plane_redshift=source_plane_redshift,
         )
-        source_triangles = triangles.with_vertices(source_plane_grid.array)
-
-        indexes = source_triangles.containing_indices(shape=shape)
-
-        return triangles.for_indexes(indexes=indexes)
+        return triangles.with_vertices(source_plane_grid.array)
 
     def steps(
         self,
@@ -278,12 +338,15 @@ def steps(
         """
         initial_triangles = self.initial_triangles
         for number in range(self.n_steps):
-            kept_triangles = self._filtered_triangles(
+            source_triangles = self._source_triangles(
                 tracer=tracer,
                 triangles=initial_triangles,
                 source_plane_redshift=source_plane_redshift,
-                shape=shape,
             )
+
+            indexes = source_triangles.containing_indices(shape=shape)
+            kept_triangles = initial_triangles.for_indexes(indexes=indexes)
+
             neighbourhood = kept_triangles
             for _ in range(self.neighbor_degree):
                 neighbourhood = neighbourhood.neighborhood()
@@ -296,6 +359,7 @@ def steps(
                 filtered_triangles=kept_triangles,
                 neighbourhood=neighbourhood,
                 up_sampled=up_sampled,
+                source_triangles=source_triangles,
             )
 
             initial_triangles = up_sampled
diff --git a/autolens/point/solver/step.py b/autolens/point/solver/step.py
index 0f06aa568..c25909677 100644
--- a/autolens/point/solver/step.py
+++ b/autolens/point/solver/step.py
@@ -5,7 +5,7 @@
 from autoarray.numpy_wrapper import register_pytree_node_class
 
 try:
-    from autoarray.structures.triangles.jax_array import ArrayTriangles
+    from autoarray.structures.triangles.array.jax_array import ArrayTriangles
 except ImportError:
     from autoarray.structures.triangles.array import ArrayTriangles
 
@@ -38,6 +38,7 @@ class Step:
     filtered_triangles: aa.AbstractTriangles
     neighbourhood: aa.AbstractTriangles
     up_sampled: aa.AbstractTriangles
+    source_triangles: aa.AbstractTriangles
 
     def tree_flatten(self):
         return (
diff --git a/autolens/point/visualise.py b/autolens/point/visualise.py
index 80de36d89..2bf253701 100644
--- a/autolens/point/visualise.py
+++ b/autolens/point/visualise.py
@@ -24,14 +24,34 @@ def visualise(step: Step):
     plt.show()
 
 
-def plot_triangles(triangles, color="black"):
+def plot_triangles(triangles, color="black", title="Triangles", point=None):
     plt.figure(figsize=(8, 8))
     for triangle in triangles:
         triangle = np.append(triangle, [triangle[0]], axis=0)
         plt.plot(triangle[:, 0], triangle[:, 1], "o-", color=color)
 
+    if point:
+        plt.plot(point[0], point[1], "x", color="red")
+
+    plt.xlabel("X")
+    plt.ylabel("Y")
+    plt.title(title)
+    plt.gca().set_aspect("equal", adjustable="box")
+    plt.show()
+
+
+def plot_triangles_compare(triangles_a, triangles_b, number=None):
+    plt.figure(figsize=(8, 8))
+    for triangle in triangles_a:
+        triangle = np.append(triangle, [triangle[0]], axis=0)
+        plt.plot(triangle[:, 0], triangle[:, 1], "o-", color="red")
+
+    for triangle in triangles_b:
+        triangle = np.append(triangle, [triangle[0]], axis=0)
+        plt.plot(triangle[:, 0], triangle[:, 1], "o-", color="blue")
+
     plt.xlabel("X")
     plt.ylabel("Y")
-    plt.title(f"Triangles")
+    plt.title("Triangles" + f" {number}" if number is not None else "")
     plt.gca().set_aspect("equal", adjustable="box")
     plt.show()
diff --git a/docs/overview/overview_1_start_here.rst b/docs/overview/overview_1_start_here.rst
index d845e8aab..1b8103a56 100644
--- a/docs/overview/overview_1_start_here.rst
+++ b/docs/overview/overview_1_start_here.rst
@@ -361,7 +361,7 @@ object.
         exposure_time=300.0,
         background_sky_level=1.0,
         psf=al.Kernel2D.from_gaussian(shape_native=(11, 11), sigma=0.1, pixel_scales=0.05),
-        add_poisson_noise=True,
+        add_poisson_noise_to_data=True,
     )
 
 Once we have a simulator, we can use it to create an imaging dataset which consists of an image, noise-map and 
diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py
index 50809b671..34917ccfc 100644
--- a/test_autolens/analysis/test_result.py
+++ b/test_autolens/analysis/test_result.py
@@ -227,7 +227,7 @@ def test__image_plane_multiple_image_positions(analysis_imaging_7x7):
 
     multiple_images = result.image_plane_multiple_image_positions
 
-    assert pytest.approx((0.968719, 0.366210), 1.0e-4) in multiple_images.in_list
+    assert pytest.approx((0.968719, 0.366210), 1.0e-2) in multiple_images.in_list
 
 
 def test__positions_threshold_from(analysis_imaging_7x7):
@@ -247,9 +247,9 @@ def test__positions_threshold_from(analysis_imaging_7x7):
 
     result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
 
-    assert result.positions_threshold_from() == pytest.approx(1.1001488121, 1.0e-4)
+    assert result.positions_threshold_from() == pytest.approx(0.930414842576, 1.0e-4)
     assert result.positions_threshold_from(factor=5.0) == pytest.approx(
-        5.5007440609, 1.0e-4
+        4.652074212, 1.0e-4
     )
     assert result.positions_threshold_from(minimum_threshold=10.0) == pytest.approx(
         10.0, 1.0e-4
@@ -291,6 +291,39 @@ def test__positions_likelihood_from(analysis_imaging_7x7):
     assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
 
 
+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.PositionsLHPenalty)
+    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
 ):
diff --git a/test_autolens/imaging/test_simulate_and_fit_imaging.py b/test_autolens/imaging/test_simulate_and_fit_imaging.py
index 8e1f2759f..06043ee84 100644
--- a/test_autolens/imaging/test_simulate_and_fit_imaging.py
+++ b/test_autolens/imaging/test_simulate_and_fit_imaging.py
@@ -25,7 +25,7 @@ def test__perfect_fit__chi_squared_0():
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.noise_map = al.Array2D.ones(
@@ -99,7 +99,7 @@ def test__perfect_fit__chi_squared_0__use_grid_iterate_to_simulate_and_fit():
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.noise_map = al.Array2D.ones(
@@ -217,7 +217,7 @@ def test__simulate_imaging_data_and_fit__linear_light_profiles_agree_with_standa
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.sub_size = 1
@@ -324,7 +324,7 @@ def test__simulate_imaging_data_and_fit__linear_light_profiles_and_pixelization(
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.sub_size = 1
@@ -477,7 +477,7 @@ def test__simulate_imaging_data_and_fit__linear_light_profiles_and_pixelization_
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.noise_map = al.Array2D.ones(
@@ -621,7 +621,7 @@ def test__simulate_imaging_data_and_fit__complex_fit_compare_mapping_matrix_w_ti
     source_1 = al.Galaxy(redshift=0.5, bulge=al.lp.Sersic(centre=(0.3, 0.3)))
     tracer = al.Tracer(galaxies=[lens_0, lens_1, lens_2, source_0, source_1])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.sub_size = 2
@@ -715,7 +715,7 @@ def test__perfect_fit__chi_squared_0__non_uniform_over_sampling():
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.noise_map = al.Array2D.ones(
@@ -803,7 +803,7 @@ def test__fit_figure_of_merit__mge_mass_model(masked_imaging_7x7, masked_imaging
     )
     tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
 
-    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise=False)
+    dataset = al.SimulatorImaging(exposure_time=300.0, psf=psf, add_poisson_noise_to_data=False)
 
     dataset = dataset.via_tracer_from(tracer=tracer, grid=grid)
     dataset.noise_map = al.Array2D.ones(
diff --git a/test_autolens/point/fit/positions/image/test_abstract.py b/test_autolens/point/fit/positions/image/test_abstract.py
index f82779813..89cd0c3df 100644
--- a/test_autolens/point/fit/positions/image/test_abstract.py
+++ b/test_autolens/point/fit/positions/image/test_abstract.py
@@ -77,8 +77,8 @@ def test__multi_plane_position_solving():
     )
 
     assert fit_0.model_data[0, 0] == pytest.approx(
-        scaling_factor * fit_1.model_data[0, 0], 1.0e-1
+        scaling_factor * fit_1.model_data[1, 0], 1.0e-1
     )
-    assert fit_0.model_data[0, 1] == pytest.approx(
+    assert fit_0.model_data[1, 1] == pytest.approx(
         scaling_factor * fit_1.model_data[0, 1], 1.0e-1
     )
diff --git a/test_autolens/point/model/test_plotter_interface_point.py b/test_autolens/point/model/test_plotter_interface_point.py
index f018daf32..7d29802e0 100644
--- a/test_autolens/point/model/test_plotter_interface_point.py
+++ b/test_autolens/point/model/test_plotter_interface_point.py
@@ -13,20 +13,16 @@ def make_plotter_interface_plotter_setup():
     return path.join("{}".format(directory), "files")
 
 
-def test__fit_point(
-    fit_point_dataset_x2_plane, include_2d_all, plot_path, plot_patch
-):
+def test__fit_point(fit_point_dataset_x2_plane, include_2d_all, plot_path, plot_patch):
     if os.path.exists(plot_path):
         shutil.rmtree(plot_path)
 
     plotter_interface = PlotterInterfacePoint(image_path=plot_path)
 
-    plotter_interface.fit_point(
-        fit=fit_point_dataset_x2_plane, during_analysis=False
-    )
+    plotter_interface.fit_point(fit=fit_point_dataset_x2_plane, during_analysis=False)
 
     assert path.join(plot_path, "subplot_fit.png") in plot_patch.paths
 
     plot_path = path.join(plot_path, "fit_dataset")
 
-    assert path.join(plot_path, "fit_point_positions.png") in plot_patch.paths
\ No newline at end of file
+    assert path.join(plot_path, "fit_point_positions.png") in plot_patch.paths
diff --git a/test_autolens/point/triangles/test_solver.py b/test_autolens/point/triangles/test_solver.py
index c75807d45..0b8625de8 100644
--- a/test_autolens/point/triangles/test_solver.py
+++ b/test_autolens/point/triangles/test_solver.py
@@ -1,9 +1,14 @@
 from typing import Tuple
 
+import numpy as np
 import pytest
 
 import autolens as al
 import autogalaxy as ag
+from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles
+from autoarray.structures.triangles.coordinate_array.jax_coordinate_array import (
+    CoordinateArrayTriangles as JAXTriangles,
+)
 from autolens.mock import NullTracer
 from autolens.point.solver import PointSolver
 
@@ -70,12 +75,39 @@ def test_trivial(
     assert coordinates[0] == pytest.approx(source_plane_coordinate, abs=1.0e-1)
 
 
-def test_real_example(grid, tracer):
-    solver = PointSolver.for_grid(
+def triangle_set(triangles):
+    return {
+        tuple(sorted([tuple(np.round(pair, 4)) for pair in triangle]))
+        for triangle in triangles.triangles.tolist()
+        if not np.isnan(triangle).any()
+    }
+
+
+def test_real_example_jax(grid, tracer):
+    jax_solver = PointSolver.for_grid(
+        grid=grid,
+        pixel_scale_precision=0.001,
+        array_triangles_cls=JAXTriangles,
+    )
+
+    result = jax_solver.solve(
+        tracer=tracer,
+        source_plane_coordinate=(0.07, 0.07),
+    )
+
+    assert len(result) == 5
+
+
+def test_real_example_normal(grid, tracer):
+    jax_solver = PointSolver.for_grid(
         grid=grid,
         pixel_scale_precision=0.001,
+        array_triangles_cls=CoordinateArrayTriangles,
     )
 
-    result = solver.solve(tracer=tracer, source_plane_coordinate=(0.07, 0.07))
+    result = jax_solver.solve(
+        tracer=tracer,
+        source_plane_coordinate=(0.07, 0.07),
+    )
 
     assert len(result) == 5
diff --git a/test_autolens/point/triangles/test_solver_jax.py b/test_autolens/point/triangles/test_solver_jax.py
index 9f5ac551f..0cbcae548 100644
--- a/test_autolens/point/triangles/test_solver_jax.py
+++ b/test_autolens/point/triangles/test_solver_jax.py
@@ -9,12 +9,13 @@
 from autolens import PointSolver, Tracer
 
 try:
-    from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles
-except ImportError:
-    from autoarray.structures.triangles.jax_coordinate_array import (
+    from autoarray.structures.triangles.coordinate_array.jax_coordinate_array import (
         CoordinateArrayTriangles,
     )
 
+except ImportError:
+    from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles
+
 from autolens.mock import NullTracer
 
 pytest.importorskip("jax")