Repository: Jammy2211/PyAutoGalaxy (language: python)

Issue:
Feature/ellipse masking
This pull request ensures that masking works for ellipse fitting correctly.

Previously, if a pixel evaluated for ellipse fitting computed one or more values for a neighboring pixel that was masked (and therefore has a value we distrust), it was still included in the likelihood calculation.

This PR ensures this behaviour no longer occurs and that any such pixel is omitted for all calculations.


Discussion / hints:

Fix patch:
diff --git a/autogalaxy/__init__.py b/autogalaxy/__init__.py
index 7f0fd842..26df5df9 100644
--- a/autogalaxy/__init__.py
+++ b/autogalaxy/__init__.py
@@ -52,8 +52,6 @@
 
 from .analysis.adapt_images.adapt_images import AdaptImages
 from .analysis.adapt_images.adapt_image_maker import AdaptImageMaker
-from .analysis.maker import FitMaker
-from .analysis.preloads import Preloads
 from . import aggregator as agg
 from . import exc
 from . import plot
@@ -116,4 +114,4 @@
 
 conf.instance.register(__file__)
 
-__version__ = "2024.5.16.0"
+__version__ = "2024.9.21.2"
diff --git a/autogalaxy/aggregator/agg_util.py b/autogalaxy/aggregator/agg_util.py
index 58b818d3..205944d4 100644
--- a/autogalaxy/aggregator/agg_util.py
+++ b/autogalaxy/aggregator/agg_util.py
@@ -68,102 +68,3 @@ def adapt_images_from(
         adapt_images_list.append(adapt_images)
 
     return adapt_images_list
-
-
-def mesh_grids_of_planes_list_from(
-    fit: af.Fit, total_fits: int, use_preloaded_grid: bool
-) -> List[Optional[aa.Grid2D]]:
-    """
-    Returns the image-plane pixelization grid(s) used by the fit.
-
-    A subset of image-mesh objects (e.g. `Hilbert`, `KMeans`) create the grid of points that act as the mesh
-    centres (e.g. the centers of Voronoi cells) in the image-plane. For lensing calculations this may then be
-    traced to the source-plane to form the pixelization.
-
-    This calculation can depend on the library used to determine the image-plane grid and may have a random element
-    associated with it. This means that performing an analysis on a super computer and then downloading the results
-    for inspection on your laptop may produce different image-plane grids, changing the result and quantities like
-    the `log_likelihood`.
-
-    By storing this grid as a result in the `files` folder and loading it via the database before creating a fit
-    this possible mismatch is removed, ensuring results on a local computer are identical to those computer elsewhere.
-
-    Parameters
-    ----------
-    fit
-        A `PyAutoFit` `Fit` object which contains the results of a model-fit as an entry in a sqlite database.
-    total_fits
-        The total number of `Analysis` objects summed to create the fit.
-    use_preloaded_grid
-        Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This grid
-        may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is used
-        as the fit.
-
-    Returns
-    -------
-    The list of image-plane mesh centres used when creating the overall fit.
-    """
-
-    if use_preloaded_grid:
-        if not fit.children:
-            return [fit.value(name="preload_mesh_grids_of_planes")]
-        else:
-            try:
-                return fit.child_values(name="preload_mesh_grids_of_planes")
-            except AttributeError:
-                return [None] * total_fits
-    else:
-        return [None] * total_fits
-
-
-def preloads_from(
-    preloads_cls,
-    use_preloaded_grid: bool,
-    mesh_grids_of_planes: List,
-    use_w_tilde: Optional[bool] = False,
-) -> aa.Preloads:
-    """
-    Returns a `Preloads` object associated with a fit loaded via the database.
-
-    When loading results via the database, the preloads class may have certain attributes associated with it
-    in order to perform the fit. The main purpose is to use the same image-plane mesh centres for pixelization where
-    the mesh is computed in the image-plane (see `agg_util.mesh_grids_of_planes_list_from`).
-
-    The preloads may also switch off `w_tilde` so fits are computed faster locally as they do not need to recompute
-    w_tilde.
-
-    Parameters
-    ----------
-    preloads_cls
-        The `Preloads` object used to create the preloads (this varies across
-        projects like `autogalaxy` and `autolens`).
-    use_preloaded_grid
-        Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This grid
-        may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is used
-        as the fit.
-    mesh_grids_of_planes
-        The list of image-plane mesh centres used when creating the overall fit which are associated with the
-        preloads.
-    use_w_tilde
-        Whether to use the w-tilde formalism when recomputing fits.
-
-    Returns
-    -------
-    The preloads object associated with the fit loaded via the database.
-    """
-    preloads = preloads_cls()
-
-    if use_preloaded_grid:
-        if mesh_grids_of_planes is not None:
-            preloads = preloads_cls(
-                image_plane_mesh_grid_pg_list=mesh_grids_of_planes,
-                use_w_tilde=use_w_tilde,
-            )
-
-            if len(preloads.image_plane_mesh_grid_pg_list) == 2:
-                if type(preloads.image_plane_mesh_grid_pg_list[1]) != list:
-                    preloads.image_plane_mesh_grid_pg_list[1] = [
-                        preloads.image_plane_mesh_grid_pg_list[1]
-                    ]
-
-    return preloads
diff --git a/autogalaxy/aggregator/ellipse/fit_ellipse.py b/autogalaxy/aggregator/ellipse/fit_ellipse.py
index 511422eb..031a7c24 100644
--- a/autogalaxy/aggregator/ellipse/fit_ellipse.py
+++ b/autogalaxy/aggregator/ellipse/fit_ellipse.py
@@ -53,15 +53,14 @@ def _fit_ellipse_from(
     try:
         multipole_list_list = _multipoles_from(fit=fit, instance=instance)
     except AttributeError:
-        multipole_list_list = [[None] * len(ellipse_list_list)] * len(
-            ellipse_list_list[0]
-        )
+        multipole_list_list = [[None for i in ellipse_list_list[0]]]
 
     fit_dataset_list = []
 
     for dataset, ellipse_list, multipole_lists in zip(
         dataset_list, ellipse_list_list, multipole_list_list
     ):
+
         for ellipse, multipole_list in zip(ellipse_list, multipole_lists):
             fit_dataset_list.append(
                 FitEllipse(
diff --git a/autogalaxy/aggregator/imaging/fit_imaging.py b/autogalaxy/aggregator/imaging/fit_imaging.py
index e805df1e..8ba62ae7 100644
--- a/autogalaxy/aggregator/imaging/fit_imaging.py
+++ b/autogalaxy/aggregator/imaging/fit_imaging.py
@@ -7,8 +7,6 @@
 import autofit as af
 import autoarray as aa
 
-from autogalaxy.analysis.preloads import Preloads
-
 from autogalaxy.aggregator.imaging.imaging import _imaging_from
 from autogalaxy.aggregator.galaxies import _galaxies_from
 from autogalaxy.aggregator.dataset_model import _dataset_model_from
@@ -19,7 +17,6 @@ def _fit_imaging_from(
     fit: af.Fit,
     instance: Optional[af.ModelInstance] = None,
     settings_inversion: aa.SettingsInversion = None,
-    use_preloaded_grid: bool = True,
 ) -> List[FitImaging]:
     """
     Returns a list of `FitImaging` objects from a `PyAutoFit` sqlite database `Fit` object.
@@ -51,10 +48,6 @@ def _fit_imaging_from(
         randomly from the PDF).
     settings_inversion
         Optionally overwrite the `SettingsInversion` of the `Inversion` object that is created from the fit.
-    use_preloaded_grid
-        Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This grid
-        may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is used
-        as the fit.
     """
 
     from autogalaxy.imaging.fit_imaging import FitImaging
@@ -69,25 +62,14 @@ def _fit_imaging_from(
 
     settings_inversion = settings_inversion or fit.value(name="settings_inversion")
 
-    mesh_grids_of_planes_list = agg_util.mesh_grids_of_planes_list_from(
-        fit=fit, total_fits=len(dataset_list), use_preloaded_grid=use_preloaded_grid
-    )
-
     fit_dataset_list = []
 
-    for dataset, galaxies, dataset_model, adapt_images, mesh_grids_of_planes in zip(
+    for dataset, galaxies, dataset_model, adapt_images in zip(
         dataset_list,
         galaxies_list,
         dataset_model_list,
         adapt_images_list,
-        mesh_grids_of_planes_list,
     ):
-        preloads = agg_util.preloads_from(
-            preloads_cls=Preloads,
-            use_preloaded_grid=use_preloaded_grid,
-            mesh_grids_of_planes=mesh_grids_of_planes,
-            use_w_tilde=False,
-        )
 
         fit_dataset_list.append(
             FitImaging(
@@ -96,7 +78,6 @@ def _fit_imaging_from(
                 dataset_model=dataset_model,
                 adapt_images=adapt_images,
                 settings_inversion=settings_inversion,
-                preloads=preloads,
             )
         )
 
@@ -108,7 +89,6 @@ def __init__(
         self,
         aggregator: af.Aggregator,
         settings_inversion: Optional[aa.SettingsInversion] = None,
-        use_preloaded_grid: bool = True,
     ):
         """
         Interfaces with an `PyAutoFit` aggregator object to create instances of `FitImaging` objects from the results
@@ -142,15 +122,10 @@ def __init__(
             A `PyAutoFit` aggregator object which can load the results of model-fits.
         settings_inversion
             Optionally overwrite the `SettingsInversion` of the `Inversion` object that is created from the fit.
-        use_preloaded_grid
-            Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This
-            grid may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is
-            used as the fit.
         """
         super().__init__(aggregator=aggregator)
 
         self.settings_inversion = settings_inversion
-        self.use_preloaded_grid = use_preloaded_grid
 
     def object_via_gen_from(
         self, fit, instance: Optional[af.ModelInstance] = None
@@ -172,5 +147,4 @@ def object_via_gen_from(
             fit=fit,
             instance=instance,
             settings_inversion=self.settings_inversion,
-            use_preloaded_grid=self.use_preloaded_grid,
         )
diff --git a/autogalaxy/aggregator/interferometer/fit_interferometer.py b/autogalaxy/aggregator/interferometer/fit_interferometer.py
index 316e3aac..5e395927 100644
--- a/autogalaxy/aggregator/interferometer/fit_interferometer.py
+++ b/autogalaxy/aggregator/interferometer/fit_interferometer.py
@@ -7,8 +7,6 @@
 import autofit as af
 import autoarray as aa
 
-from autogalaxy.analysis.preloads import Preloads
-
 from autogalaxy.aggregator import agg_util
 from autogalaxy.aggregator.interferometer.interferometer import _interferometer_from
 from autogalaxy.aggregator.dataset_model import _dataset_model_from
@@ -20,7 +18,6 @@ def _fit_interferometer_from(
     instance: Optional[af.ModelInstance] = None,
     real_space_mask: Optional[aa.Mask2D] = None,
     settings_inversion: aa.SettingsInversion = None,
-    use_preloaded_grid: bool = True,
 ) -> List[FitInterferometer]:
     """
     Returns a list of `FitInterferometer` objects from a `PyAutoFit` sqlite database `Fit` object.
@@ -53,10 +50,6 @@ def _fit_interferometer_from(
         randomly from the PDF).
     settings_inversion
         Optionally overwrite the `SettingsInversion` of the `Inversion` object that is created from the fit.
-    use_preloaded_grid
-        Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This grid
-        may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is used
-        as the fit.
     """
     from autogalaxy.interferometer.fit_interferometer import FitInterferometer
 
@@ -73,25 +66,14 @@ def _fit_interferometer_from(
 
     settings_inversion = settings_inversion or fit.value(name="settings_inversion")
 
-    mesh_grids_of_planes_list = agg_util.mesh_grids_of_planes_list_from(
-        fit=fit, total_fits=len(dataset_list), use_preloaded_grid=use_preloaded_grid
-    )
-
     fit_dataset_list = []
 
-    for dataset, galaxies, dataset_model, adapt_images, mesh_grids_of_planes in zip(
+    for dataset, galaxies, dataset_model, adapt_images in zip(
         dataset_list,
         galaxies_list,
         dataset_model_list,
         adapt_images_list,
-        mesh_grids_of_planes_list,
     ):
-        preloads = agg_util.preloads_from(
-            preloads_cls=Preloads,
-            use_preloaded_grid=use_preloaded_grid,
-            mesh_grids_of_planes=mesh_grids_of_planes,
-            use_w_tilde=False,
-        )
 
         fit_dataset_list.append(
             FitInterferometer(
@@ -100,7 +82,6 @@ def _fit_interferometer_from(
                 dataset_model=dataset_model,
                 adapt_images=adapt_images,
                 settings_inversion=settings_inversion,
-                preloads=preloads,
             )
         )
 
@@ -112,7 +93,6 @@ def __init__(
         self,
         aggregator: af.Aggregator,
         settings_inversion: Optional[aa.SettingsInversion] = None,
-        use_preloaded_grid: bool = True,
         real_space_mask: Optional[aa.Mask2D] = None,
     ):
         """
@@ -147,15 +127,10 @@ def __init__(
             A `PyAutoFit` aggregator object which can load the results of model-fits.
         settings_inversion
             Optionally overwrite the `SettingsInversion` of the `Inversion` object that is created from the fit.
-        use_preloaded_grid
-            Certain pixelization's construct their mesh in the source-plane from a stochastic KMeans algorithm. This
-            grid may be output to hard-disk after the model-fit and loaded via the database to ensure the same grid is
-            used as the fit.
         """
         super().__init__(aggregator=aggregator)
 
         self.settings_inversion = settings_inversion
-        self.use_preloaded_grid = use_preloaded_grid
         self.real_space_mask = real_space_mask
 
     def object_via_gen_from(
@@ -178,5 +153,4 @@ def object_via_gen_from(
             fit=fit,
             instance=instance,
             settings_inversion=self.settings_inversion,
-            use_preloaded_grid=self.use_preloaded_grid,
         )
diff --git a/autogalaxy/analysis/analysis/analysis.py b/autogalaxy/analysis/analysis/analysis.py
index f084d1f1..15645c66 100644
--- a/autogalaxy/analysis/analysis/analysis.py
+++ b/autogalaxy/analysis/analysis/analysis.py
@@ -138,7 +138,7 @@ def profile_log_likelihood_function(
         An `info_dict` is also created which stores information on aspects of the model and dataset that dictate
         run times, so the profiled times can be interpreted with this context.
 
-        The results of this profiling are then output to hard-disk in the `preloads` folder of the model-fit results,
+        The results of this profiling are then output to hard-disk in the `prefiling` folder of the model-fit results,
         which they can be inspected to ensure run-times are as expected.
 
         Parameters
@@ -191,33 +191,33 @@ def profile_log_likelihood_function(
 
         try:
             info_dict["image_pixels"] = self.dataset.grids.uniform.shape_slim
-            info_dict[
-                "sub_total_light_profiles"
-            ] = self.dataset.grids.uniform.over_sampler.sub_total
+            info_dict["sub_total_light_profiles"] = (
+                self.dataset.grids.uniform.over_sampler.sub_total
+            )
         except AttributeError:
             pass
 
         if fit.model_obj.has(cls=aa.Pixelization):
             info_dict["use_w_tilde"] = fit.inversion.settings.use_w_tilde
             try:
-                info_dict[
-                    "sub_total_pixelization"
-                ] = self.dataset.grids.pixelization.over_sampler.sub_total
+                info_dict["sub_total_pixelization"] = (
+                    self.dataset.grids.pixelization.over_sampler.sub_total
+                )
             except AttributeError:
                 pass
-            info_dict[
-                "use_positive_only_solver"
-            ] = fit.inversion.settings.use_positive_only_solver
-            info_dict[
-                "force_edge_pixels_to_zeros"
-            ] = fit.inversion.settings.force_edge_pixels_to_zeros
+            info_dict["use_positive_only_solver"] = (
+                fit.inversion.settings.use_positive_only_solver
+            )
+            info_dict["force_edge_pixels_to_zeros"] = (
+                fit.inversion.settings.force_edge_pixels_to_zeros
+            )
             info_dict["use_w_tilde_numpy"] = fit.inversion.settings.use_w_tilde_numpy
             info_dict["source_pixels"] = len(fit.inversion.reconstruction)
 
             if hasattr(fit.inversion, "w_tilde"):
-                info_dict[
-                    "w_tilde_curvature_preload_size"
-                ] = fit.inversion.w_tilde.curvature_preload.shape[0]
+                info_dict["w_tilde_curvature_preload_size"] = (
+                    fit.inversion.w_tilde.curvature_preload.shape[0]
+                )
 
         self.output_profiling_info(
             paths=paths, run_time_dict=run_time_dict, info_dict=info_dict
diff --git a/autogalaxy/analysis/analysis/dataset.py b/autogalaxy/analysis/analysis/dataset.py
index da8a3571..ac5e957e 100644
--- a/autogalaxy/analysis/analysis/dataset.py
+++ b/autogalaxy/analysis/analysis/dataset.py
@@ -1,17 +1,13 @@
 import copy
 import logging
 from typing import Optional, Union
-import os
 
-from autoconf import conf
 from autoconf.dictable import to_dict, output_to_json
 import autofit as af
 import autoarray as aa
 
 from autogalaxy.analysis.adapt_images.adapt_image_maker import AdaptImageMaker
 from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages
-from autogalaxy.analysis.maker import FitMaker
-from autogalaxy.analysis.preloads import Preloads
 from autogalaxy.cosmology.lensing import LensingCosmology
 from autogalaxy.cosmology.wrap import Planck15
 from autogalaxy.analysis.analysis.analysis import Analysis
@@ -62,18 +58,8 @@ def __init__(
 
         self.settings_inversion = settings_inversion or aa.SettingsInversion()
 
-        self.preloads = self.preloads_cls()
-
         self.title_prefix = title_prefix
 
-    @property
-    def preloads_cls(self):
-        return Preloads
-
-    @property
-    def fit_maker_cls(self):
-        return FitMaker
-
     @property
     def adapt_images(self):
         if self._adapt_images is not None:
@@ -82,67 +68,6 @@ def adapt_images(self):
         if self.adapt_image_maker is not None:
             return self.adapt_image_maker.adapt_images
 
-    def set_preloads(self, paths: af.DirectoryPaths, model: af.Collection):
-        """
-        It is common for the model to have components whose parameters are all fixed, and thus the way that component
-        fits the data does not change. For example, if all parameter associated with the light profiles of galaxies
-        in the model are fixed, the image generated from these galaxies will not change irrespective of the model
-        parameters chosen by the non-linear search.
-
-        Preloading exploits this to speed up the log likelihood function, by inspecting the model and storing in memory
-        quantities that do not change. For the example above, the image of all galaxies would be stored in memory and
-        to perform every fit in the `log_likelihood_funtion`.
-
-        This function sets up all preload quantities, which are described fully in the `preloads` modules. This
-        occurs directly before the non-linear search begins, to ensure the model parameterization is fixed.
-
-        Parameters
-        ----------
-        paths
-            The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored,
-            visualization and the pickled objects used by the aggregator output by this function.
-        model
-            The PyAutoFit model object, which includes model components representing the galaxies that are fitted to
-            the imaging data.
-        """
-
-        logger.info(
-            "PRELOADS - Setting up preloads, may take a few minutes for fits using an inversion."
-        )
-
-        self.preloads = self.preloads_cls()
-
-        settings_inversion_original = copy.copy(self.settings_inversion)
-
-        self.settings_inversion.image_mesh_min_mesh_pixels_per_pixel = None
-        self.settings_inversion.image_mesh_adapt_background_percent_threshold = None
-
-        fit_maker = self.fit_maker_cls(model=model, fit_from=self.fit_from)
-
-        fit_0 = fit_maker.fit_via_model_from(unit_value=0.45)
-        fit_1 = fit_maker.fit_via_model_from(unit_value=0.55)
-
-        if fit_0 is None or fit_1 is None:
-            self.preloads = self.preloads_cls(failed=True)
-
-            self.settings_inversion = settings_inversion_original
-
-        else:
-            self.preloads = self.preloads_cls.setup_all_via_fits(
-                fit_0=fit_0, fit_1=fit_1
-            )
-
-            if conf.instance["general"]["test"]["check_preloads"]:
-                self.preloads.check_via_fit(fit=fit_0)
-
-            self.settings_inversion = settings_inversion_original
-
-        if isinstance(paths, af.DatabasePaths):
-            return
-
-        os.makedirs(paths.profile_path, exist_ok=True)
-        self.preloads.output_info_to_summary(file_path=paths.profile_path)
-
     def save_attributes(self, paths: af.DirectoryPaths):
         """
         Before the model-fit via the non-linear search begins, this routine saves attributes of the `Analysis` object
diff --git a/autogalaxy/analysis/maker.py b/autogalaxy/analysis/maker.py
deleted file mode 100644
index 154042ee..00000000
--- a/autogalaxy/analysis/maker.py
+++ /dev/null
@@ -1,128 +0,0 @@
-import logging
-from typing import Callable, Union
-
-import autoarray as aa
-import autofit as af
-from autoconf import conf
-from autofit.exc import PriorLimitException
-from autogalaxy.analysis.preloads import Preloads
-
-logger = logging.getLogger(__name__)
-
-logger.setLevel(level="INFO")
-
-
-class FitMaker:
-    def __init__(self, model: af.Collection, fit_from: Callable):
-        """
-        Makes fits using an input PyAutoFit `model`, where the parameters of the model are drawn from its prior. This
-        uses an input `fit_from`, which given an `instance` of the model creates the fit object.
-
-        This is used for implicit preloading in the `Analysis` classes, whereby the created fits are compared against
-        one another to determine whether certain components of the analysis can be preloaded.
-
-        This includes functionality for creating the fit via the model in different ways, so that if certain
-        models are ill-defined another is used instead.
-
-        Parameters
-        ----------
-        model
-            A **PyAutoFit** model object which via its parameters and their priors can created instances of the model.
-        fit_from
-            A function which given the instance of the model creates a `Fit` object.
-        """
-
-        self.model = model
-        self.fit_from = fit_from
-
-    @property
-    def preloads_cls(self):
-        return Preloads
-
-    def fit_via_model_from(
-        self, unit_value: float
-    ) -> Union[aa.FitImaging, aa.FitInterferometer]:
-        """
-        Create a fit via the model.
-
-        This first tries to compute the fit from the input `unit_value`, where the `unit_value` defines unit hyper
-        cube values of each parameter's prior in the model, used to map each value to physical values for the fit.
-
-        If this model fit produces an `Exception` because the parameter combination does not fit the data accurately,
-        a sequence of random fits are instead used into an exception is not returned. However, if the number
-        of `preload_attempts` defined in the configuration files is exceeded a None is returned.
-
-        Parameters
-        ----------
-        unit_value
-            The unit hyper cube values of each parameter's prior in the model, used to map each value to physical
-            values for the fit.
-
-        Returns
-        -------
-        fit
-            A fit object where an instance of the model has been fitted to the data.
-        """
-        try:
-            try:
-                return self.fit_unit_instance_from(unit_value=unit_value)
-            except IndexError as e:
-                raise Exception from e
-        except (Exception, PriorLimitException):
-            return self.fit_random_instance_from()
-
-    def fit_unit_instance_from(
-        self, unit_value: float
-    ) -> Union[aa.FitImaging, aa.FitInterferometer]:
-        """
-        Create a fit via the model using an input `unit_value`, where the `unit_value` defines unit hyper
-        cube values of each parameter's prior in the model, used to map each value to physical values for the fit.
-
-        Parameters
-        ----------
-        unit_value
-            The unit hyper cube values of each parameter's prior in the model, used to map each value to physical
-            values for the fit.
-
-        Returns
-        -------
-        fit
-            A fit object where an instance of the model has been fitted to the data.
-        """
-        instance = self.model.instance_from_unit_vector(
-            unit_vector=[unit_value] * self.model.prior_count, ignore_prior_limits=True
-        )
-
-        fit = self.fit_from(
-            instance=instance,
-        )
-        fit.figure_of_merit
-        return fit
-
-    def fit_random_instance_from(self) -> Union[aa.FitImaging, aa.FitInterferometer]:
-        """
-        Create a fit via the model by guessing a sequence of random fits until an exception is not returned. If
-        the number of `preload_attempts` defined in the configuration files is exceeded a None is returned.
-
-        Returns
-        -------
-        fit
-            A fit object where an instance of the model has been fitted to the data.
-        """
-
-        preload_attempts = conf.instance["general"]["analysis"]["preload_attempts"]
-
-        for i in range(preload_attempts):
-            try:
-                instance = self.model.random_instance(ignore_prior_limits=True)
-
-                fit = self.fit_from(
-                    instance=instance,
-                )
-
-                fit.figure_of_merit
-
-                return fit
-
-            except Exception as e:
-                continue
diff --git a/autogalaxy/analysis/preloads.py b/autogalaxy/analysis/preloads.py
deleted file mode 100644
index b1d513e8..00000000
--- a/autogalaxy/analysis/preloads.py
+++ /dev/null
@@ -1,213 +0,0 @@
-import logging
-import numpy as np
-from os import path
-from typing import Dict, Optional, List
-
-import autofit as af
-import autoarray as aa
-
-logger = logging.getLogger(__name__)
-
-logger.setLevel(level="INFO")
-
-
-class Preloads(aa.Preloads):
-    def __init__(
-        self,
-        w_tilde: Optional[aa.WTildeImaging] = None,
-        use_w_tilde: Optional[bool] = None,
-        blurred_image: Optional[aa.Array2D] = None,
-        image_plane_mesh_grid_pg_list: Optional[List[List[aa.Grid2D]]] = None,
-        relocated_grid: Optional[aa.Grid2D] = None,
-        mapper_list: Optional[aa.AbstractMapper] = None,
-        mapper_galaxy_dict: Optional[Dict[aa.AbstractMapper, "Galaxy"]] = None,
-        linear_func_operated_mapping_matrix_dict=None,
-        data_linear_func_matrix_dict=None,
-        mapper_operated_mapping_matrix_dict=None,
-        operated_mapping_matrix: Optional[np.ndarray] = None,
-        curvature_matrix: Optional[np.ndarray] = None,
-        regularization_matrix: Optional[np.ndarray] = None,
-        log_det_regularization_matrix_term: Optional[float] = None,
-        traced_mesh_grids_list_of_planes=None,
-        image_plane_mesh_grid_list=None,
-        failed=False,
-    ):
-        """
-        Class which offers a concise API for settings up the preloads, which before a model-fit are set up via
-        a comparison of two fits using two different models. If a quantity in these two fits is identical, it does
-        not change thoughout the model-fit and can therefore be preloaded to avoid computation, speeding up
-        the analysis.
-
-        For example, the image-plane source-plane pixelization grid (which may be computationally expensive to compute
-        via a KMeans algorithm) does not change for the majority of model-fits, because the associated model parameters
-        are fixed. Preloading avoids rerruning the KMeans algorithm for every model fitted, by preloading it in memory
-        and using this preload in every fit.
-
-        Parameters
-        ----------
-        blurred_image
-            The preloaded array of values containing the blurred image of a model fit (e.g. that light profile of
-            every galaxy in the model). This can be preloaded when no light profiles in the model vary.
-        w_tilde
-            A class containing values that enable an inversion's linear algebra to use the w-tilde formalism. This can
-            be preloaded when no component of the model changes the noise map (e.g. galaxies are fixed).
-        use_w_tilde
-            Whether to use the w tilde formalism, which superseeds the value in `SettingsInversions` such that w tilde
-            will be disabled for model-fits it is not applicable (e.g. because the noise-map changes).
-        traced_grids_of_planes_for_inversion
-            The two dimensional grids corresponding to the traced grids in a lens fit. This can be preloaded when no
-             mass profiles in the model vary.
-        image_plane_mesh_grid_pg_list
-            The two dimensional grids corresponding to the sparse image plane grids in a lens fit, that is ray-traced to
-            the source plane to form the source pixelization. This can be preloaded when no pixelizations in the model
-            vary.
-        relocated_grid
-            The two dimensional grids corresponding to the grid that has had its border pixels relocated for a
-            pixelization in a lens fit. This can be preloaded when no mass profiles in the model vary.
-        mapper_list
-            The mapper of a fit, which preloading avoids recalculation of the mapping matrix and image to source
-            pixel mappings. This can be preloaded when no pixelizations in the model vary.
-        operated_mapping_matrix
-            A matrix containing the mappings between PSF blurred image pixels and source pixels used in the linear
-            algebra of an inversion. This can be preloaded when no mass profiles and pixelizations in the model vary.
-
-        Returns
-        -------
-        The preloads object used to skip certain calculations in the log likelihood function.
-        """
-        super().__init__(
-            w_tilde=w_tilde,
-            use_w_tilde=use_w_tilde,
-            relocated_grid=relocated_grid,
-            image_plane_mesh_grid_pg_list=image_plane_mesh_grid_pg_list,
-            mapper_list=mapper_list,
-            linear_func_operated_mapping_matrix_dict=linear_func_operated_mapping_matrix_dict,
-            data_linear_func_matrix_dict=data_linear_func_matrix_dict,
-            mapper_operated_mapping_matrix_dict=mapper_operated_mapping_matrix_dict,
-            operated_mapping_matrix=operated_mapping_matrix,
-            curvature_matrix=curvature_matrix,
-            regularization_matrix=regularization_matrix,
-            log_det_regularization_matrix_term=log_det_regularization_matrix_term,
-            traced_mesh_grids_list_of_planes=traced_mesh_grids_list_of_planes,
-            image_plane_mesh_grid_list=image_plane_mesh_grid_list,
-        )
-
-        self.mapper_galaxy_dict = mapper_galaxy_dict
-        self.blurred_image = blurred_image
-        self.failed = failed
-
-    @classmethod
-    def setup_all_via_fits(cls, fit_0, fit_1) -> "Preloads":
-        """
-        Setup the Preloads from two fits which use two different lens model of a model-fit.
-
-        Parameters
-        ----------
-        fit_0
-            The first fit corresponding to a model with a specific set of unit-values.
-        fit_1
-            The second fit corresponding to a model with a different set of unit-values.
-
-        Returns
-        -------
-        Preloads
-            Preloads which are set up based on the fit's passed in specific to a lens model.
-
-        """
-
-        preloads = cls()
-
-        if isinstance(fit_0, aa.FitImaging):
-            preloads.set_w_tilde_imaging(fit_0=fit_0, fit_1=fit_1)
-            preloads.set_blurred_image(fit_0=fit_0, fit_1=fit_1)
-
-        preloads.set_mapper_list(fit_0=fit_0, fit_1=fit_1)
-
-        if preloads.mapper_list is not None:
-            preloads.mapper_galaxy_dict = fit_0.galaxies_to_inversion.mapper_galaxy_dict
-
-        preloads.set_operated_mapping_matrix_with_preloads(fit_0=fit_0, fit_1=fit_1)
-        preloads.set_linear_func_inversion_dicts(fit_0=fit_0, fit_1=fit_1)
-        preloads.set_curvature_matrix(fit_0=fit_0, fit_1=fit_1)
-        preloads.set_regularization_matrix_and_term(fit_0=fit_0, fit_1=fit_1)
-
-        return preloads
-
-    def set_blurred_image(self, fit_0, fit_1):
-        """
-        If the `LightProfile`'s in a model are all fixed parameters their corresponding image and therefore PSF blurred
-        image do not change during the model fit and can therefore be preloaded.
-
-        This function compares the blurred image of two fit's corresponding to two model instances, and preloads
-        the blurred image if the blurred image of both fits are the same.
-
-        The preload is typically used though out search chaining pipelines, as it is common to fix the lens light for
-        the majority of model-fits.
-
-        Parameters
-        ----------
-        fit_0
-            The first fit corresponding to a model with a specific set of unit-values.
-        fit_1
-            The second fit corresponding to a model with a different set of unit-values.
-        """
-        self.blurred_image = None
-
-        precision = 1e-8
-
-        if (np.max(abs(fit_0.blurred_image - fit_1.blurred_image)) < precision) and (
-            np.sum(fit_0.blurred_image) > precision
-        ):
-            self.blurred_image = fit_0.blurred_image
-
-            logger.info(
-                "PRELOADS - Blurred image (e.g. the image of all light profiles) is preloaded for this model-fit."
-            )
-
-    def output_info_to_summary(self, file_path):
-        file_preloads = path.join(file_path, "preloads.summary")
-
-        if self.failed:
-            logger.info(
-                "PRELOADS - Preloading failed as models gave too many exceptions, preloading therefore "
-                "not used."
-            )
-
-            af.formatter.output_list_of_strings_to_file(
-                file=file_preloads, list_of_strings=["FAILED"]
-            )
-
-        else:
-            af.formatter.output_list_of_strings_to_file(
-                file=file_preloads, list_of_strings=self.info
-            )
-
-    @property
-    def info(self) -> List[str]:
-        """
-        The information on what has or has not been preloaded, which is written to the file `preloads.summary`.
-
-        Returns
-        -------
-            A list of strings containing True and False values as to whether a quantity has been preloaded.
-        """
-        line = [f"W Tilde = {self.w_tilde is not None}\n"]
-        line += [f"Use W Tilde = {self.use_w_tilde}\n\n"]
-        line += [f"Blurred Image = {np.count_nonzero(self.blurred_image) != 0}\n"]
-        line += [f"Mapper = {self.mapper_list is not None}\n"]
-        line += [
-            f"Blurred Mapping Matrix = {self.operated_mapping_matrix is not None}\n"
-        ]
-        line += [
-            f"Inversion Linear Func (Linear Light Profile) Dicts = {self.linear_func_operated_mapping_matrix_dict is not None}\n"
-        ]
-        line += [f"Curvature Matrix = {self.curvature_matrix is not None}\n"]
-        line += [
-            f"Curvature Matrix Mapper Diag = {self.curvature_matrix_mapper_diag is not None}\n"
-        ]
-        line += [f"Regularization Matrix = {self.regularization_matrix is not None}\n"]
-        line += [
-            f"Log Det Regularization Matrix Term = {self.log_det_regularization_matrix_term is not None}\n"
-        ]
-
-        return line
diff --git a/autogalaxy/config/general.yaml b/autogalaxy/config/general.yaml
index 61f8b09a..361abe73 100644
--- a/autogalaxy/config/general.yaml
+++ b/autogalaxy/config/general.yaml
@@ -1,5 +1,3 @@
-analysis:
-  preload_attempts: 250
 fits:
   flip_for_ds9: true
 grid:
@@ -11,6 +9,4 @@ inversion:
   use_border_relocator: true          # If True, by default a pixelization's border is used to relocate all pixels outside its border to the border.
 test:
   check_likelihood_function: true   # if True, when a search is resumed the likelihood of a previous sample is recalculated to ensure it is consistent with the previous run.
-  check_preloads: false
-  exception_override: false
-  preloads_check_threshold: 1.0     # If the figure of merit of a fit with and without preloads is greater than this threshold, the check preload test fails and an exception raised for a model-fit. 
+  exception_override: false
\ No newline at end of file
diff --git a/autogalaxy/ellipse/dataset_interp.py b/autogalaxy/ellipse/dataset_interp.py
index 997951af..ea165c5b 100644
--- a/autogalaxy/ellipse/dataset_interp.py
+++ b/autogalaxy/ellipse/dataset_interp.py
@@ -44,6 +44,19 @@ def points_interp(self) -> Tuple[np.ndarray, np.ndarray]:
 
         return (x, y)
 
+    @cached_property
+    def mask_interp(self) -> interpolate.RegularGridInterpolator:
+        """
+        Returns a 2D interpolation of the mask, which is used to determine whether inteprolated values use a masked
+        pixel for the interpolation and thus should not be included in a fit.
+        """
+        return interpolate.RegularGridInterpolator(
+            points=self.points_interp,
+            values=np.float64(self.dataset.data.mask),
+            bounds_error=False,
+            fill_value=0.0,
+        )
+
     @cached_property
     def data_interp(self) -> interpolate.RegularGridInterpolator:
         """
diff --git a/autogalaxy/ellipse/ellipse/ellipse.py b/autogalaxy/ellipse/ellipse/ellipse.py
index d8bf9a27..b6d1acd0 100644
--- a/autogalaxy/ellipse/ellipse/ellipse.py
+++ b/autogalaxy/ellipse/ellipse/ellipse.py
@@ -115,6 +115,9 @@ def angles_from_x0_from(self, pixel_scale: float) -> np.ndarray:
 
         They are subtracted by the `angle` of the ellipse to give the angles from the major-axis of the ellipse.
 
+        The final angle, which is 2.0 * np.pi radians, is a repeat of the first angle of zero radians, therefore
+        the final angle is removed.
+
         The number of angles computed is the minimum of 500 and the integer of the circular radius of the ellipse.
         This value is chosen to ensure that the number of angles computed matches the number of pixels in the
         data that the ellipse is being fitted to.
@@ -130,7 +133,7 @@ def angles_from_x0_from(self, pixel_scale: float) -> np.ndarray:
         """
         total_points = self.total_points_from(pixel_scale)
 
-        return np.linspace(0.0, 2.0 * np.pi, total_points)
+        return np.linspace(0.0, 2.0 * np.pi, total_points)[:-1]
 
     def ellipse_radii_from_major_axis_from(self, pixel_scale: float) -> np.ndarray:
         """
diff --git a/autogalaxy/ellipse/fit_ellipse.py b/autogalaxy/ellipse/fit_ellipse.py
index 446eeba6..6ea2ba94 100644
--- a/autogalaxy/ellipse/fit_ellipse.py
+++ b/autogalaxy/ellipse/fit_ellipse.py
@@ -80,6 +80,52 @@ def _points_from_major_axis(self) -> np.ndarray:
         """
         return self.points_from_major_axis_from()
 
+    @property
+    def mask_interp(self) -> np.ndarray:
+        """
+        Returns the mask values of the dataset that the ellipse fits, which are computed by overlaying the ellipse over
+        the 2D data and performing a 2D interpolation at discrete (y,x) coordinates on the ellipse on the dataset's
+        mask.
+
+        When an input (y,x) coordinate intepolates only unmasked values (`data.mask=False`) the intepolatred value
+        is 0.0, where if it interpolates one or a masked value (`data.mask=True`), the interpolated value is positive.
+        To mask all values which interpolate a masked value, all interpolated values above 1 and converted to `True`.
+
+        This mask is used to remove these pixels from a fit and evaluate how many ellipse points are used for each
+        ellipse fit.
+
+        The (y,x) coordinates on the ellipse where the interpolation occurs are computed in the
+        `points_from_major_axis` property of the `Ellipse` class, with the documentation describing how these points
+        are computed.
+
+        Returns
+        -------
+        The data values of the ellipse fits, computed via a 2D interpolation of where the ellipse
+        overlaps the data.
+        """
+        return self.interp.mask_interp(self._points_from_major_axis) > 0.0
+
+    @property
+    def total_points_interp(self) -> int:
+        """
+        Returns the total number of points used to interpolate the data and noise-map values of the ellipse.
+
+        For example, if the ellipse spans 10 pixels, the total number of points will be 10.
+
+        The calculation removes points if one or more interpolated value uses a masked value, meaning the interpolation
+        is not reliable.
+
+        The (y,x) coordinates on the ellipse where the interpolation occurs are computed in the
+        `points_from_major_axis` property of the `Ellipse` class, with the documentation describing how these points
+        are computed.
+
+        Returns
+        -------
+        The noise-map values of the ellipse fits, computed via a 2D interpolation of where the ellipse
+        overlaps the noise-map.
+        """
+        return self.data_interp[np.invert(self.mask_interp)].shape[0]
+
     @property
     def data_interp(self) -> aa.ArrayIrregular:
         """
@@ -90,15 +136,19 @@ def data_interp(self) -> aa.ArrayIrregular:
         `points_from_major_axis` property of the `Ellipse` class, with the documentation describing how these points
         are computed.
 
+        If the interpolation of an ellipse point uses one or more masked values, this point is not reliable, therefore
+        the value is converted to `np.nan` and not used by other fitting quantities.
+
         Returns
         -------
         The data values of the ellipse fits, computed via a 2D interpolation of where the ellipse
         overlaps the data.
         """
+        data = self.interp.data_interp(self._points_from_major_axis)
 
-        return aa.ArrayIrregular(
-            values=self.interp.data_interp(self._points_from_major_axis)
-        )
+        data[self.mask_interp] = np.nan
+
+        return aa.ArrayIrregular(values=data)
 
     @property
     def noise_map_interp(self) -> aa.ArrayIrregular:
@@ -110,14 +160,19 @@ def noise_map_interp(self) -> aa.ArrayIrregular:
         `points_from_major_axis` property of the `Ellipse` class, with the documentation describing how these points
         are computed.
 
+        If the interpolation of an ellipse point uses one or more masked values, this point is not reliable, therefore
+        the value is converted to `np.nan` and not used by other fitting quantities.
+
         Returns
         -------
         The noise-map values of the ellipse fits, computed via a 2D interpolation of where the ellipse
         overlaps the noise-map.
         """
-        return aa.ArrayIrregular(
-            values=self.interp.noise_map_interp(self._points_from_major_axis)
-        )
+        noise_map = self.interp.noise_map_interp(self._points_from_major_axis)
+
+        noise_map[self.mask_interp] = np.nan
+
+        return aa.ArrayIrregular(values=noise_map)
 
     @property
     def signal_to_noise_map_interp(self) -> aa.ArrayIrregular:
@@ -153,7 +208,7 @@ def model_data(self) -> aa.ArrayIrregular:
         return self.data_interp
 
     @property
-    def residual_map(self):
+    def residual_map(self) -> aa.ArrayIrregular:
         """
         Returns the residual-map of the fit, which is the data minus the model data and therefore the same
         as the model data.
@@ -162,6 +217,10 @@ def residual_map(self):
         -------
         The residual-map of the fit, which is the data minus the model data and therefore the same as the model data.
         """
+
+        if not self.model_data:
+            return aa.ArrayIrregular(values=np.zeros(self.model_data.shape))
+
         return aa.ArrayIrregular(values=self.model_data - np.nanmean(self.model_data))
 
     @property
@@ -177,16 +236,7 @@ def normalized_residual_map(self) -> aa.ArrayIrregular:
         -------
         The normalized residual-map of the fit, which is the residual-map divided by the noise-map.
         """
-
-        normalized_residual_map = (self.residual_map) / self.noise_map_interp
-
-        # NOTE:
-        idx = np.logical_or(
-            np.isnan(normalized_residual_map), np.isinf(normalized_residual_map)
-        )
-        normalized_residual_map[idx] = 0.0
-
-        return aa.ArrayIrregular(values=normalized_residual_map)
+        return aa.ArrayIrregular(values=self.residual_map / self.noise_map_interp)
 
     @property
     def chi_squared_map(self) -> aa.ArrayIrregular:
@@ -216,7 +266,7 @@ def chi_squared(self) -> float:
         -------
         The chi-squared of the fit.
         """
-        return float(np.sum(self.chi_squared_map))
+        return float(np.nansum(self.chi_squared_map))
 
     @property
     def noise_normalization(self):
@@ -227,7 +277,7 @@ def noise_normalization(self):
         -------
         The noise normalization term of the log likelihood.
         """
-        return np.sum(np.log(2 * np.pi * self.noise_map_interp**2.0))
+        return np.nansum(np.log(2 * np.pi * self.noise_map_interp**2.0))
 
     @property
     def log_likelihood(self):
diff --git a/autogalaxy/ellipse/model/analysis.py b/autogalaxy/ellipse/model/analysis.py
index cd51e9a5..07f92f41 100644
--- a/autogalaxy/ellipse/model/analysis.py
+++ b/autogalaxy/ellipse/model/analysis.py
@@ -217,7 +217,7 @@ def profile_log_likelihood_function(
         An `info_dict` is also created which stores information on aspects of the model and dataset that dictate
         run times, so the profiled times can be interpreted with this context.
 
-        The results of this profiling are then output to hard-disk in the `preloads` folder of the model-fit results,
+        The results of this profiling are then output to hard-disk in the `profiling` folder of the model-fit results,
         which they can be inspected to ensure run-times are as expected.
 
         Parameters
diff --git a/autogalaxy/galaxy/to_inversion.py b/autogalaxy/galaxy/to_inversion.py
index 669d62f9..76d233ed 100644
--- a/autogalaxy/galaxy/to_inversion.py
+++ b/autogalaxy/galaxy/to_inversion.py
@@ -15,7 +15,6 @@
 from autogalaxy.profiles.light.linear import LightProfileLinear
 from autogalaxy.galaxy.galaxy import Galaxy
 from autogalaxy.galaxy.galaxies import Galaxies
-from autogalaxy.analysis.preloads import Preloads
 
 
 class AbstractToInversion:
@@ -24,7 +23,6 @@ def __init__(
         dataset: Optional[Union[aa.Imaging, aa.Interferometer, aa.DatasetInterface]],
         adapt_images: Optional[AdaptImages] = None,
         settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads=Preloads(),
         run_time_dict: Optional[Dict] = None,
     ):
         """
@@ -59,8 +57,6 @@ def __init__(
             the pixelization's pixels to the brightest regions of the image.
         settings_inversion
             The settings of the inversion, which controls how the linear algebra calculation is performed.
-        preloads
-            Preloads of the inversion, which are used to speed up the linear algebra calculation.
         run_time_dict
             A dictionary of run-time values used to compute the inversion, for example the noise-map normalization.
         """
@@ -82,7 +78,6 @@ def __init__(
 
         self.settings_inversion = settings_inversion
 
-        self.preloads = preloads
         self.run_time_dict = run_time_dict
 
     @property
@@ -198,7 +193,6 @@ def __init__(
         galaxies: List[Galaxy],
         adapt_images: Optional[AdaptImages] = None,
         settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads=aa.Preloads(),
         run_time_dict: Optional[Dict] = None,
     ):
         """
@@ -234,8 +228,6 @@ def __init__(
             the pixelization's pixels to the brightest regions of the image.
         settings_inversion
             The settings of the inversion, which controls how the linear algebra calculation is performed.
-        preloads
-            Preloads of the inversion, which are used to speed up the linear algebra calculation.
         run_time_dict
             A dictionary of run-time values used to compute the inversion, for example the noise-map normalization.
         """
@@ -245,7 +237,6 @@ def __init__(
             dataset=dataset,
             adapt_images=adapt_images,
             settings_inversion=settings_inversion,
-            preloads=preloads,
             run_time_dict=run_time_dict,
         )
 
@@ -500,7 +491,6 @@ def mapper_from(
             source_plane_mesh_grid=source_plane_mesh_grid,
             image_plane_mesh_grid=image_plane_mesh_grid,
             adapt_data=adapt_galaxy_image,
-            preloads=self.preloads,
             run_time_dict=self.run_time_dict,
         )
 
@@ -592,7 +582,6 @@ def inversion(self) -> aa.AbstractInversion:
             dataset=self.dataset,
             linear_obj_list=self.linear_obj_list,
             settings=self.settings_inversion,
-            preloads=self.preloads,
             run_time_dict=self.run_time_dict,
         )
 
diff --git a/autogalaxy/imaging/fit_imaging.py b/autogalaxy/imaging/fit_imaging.py
index ff72fe5e..06f33169 100644
--- a/autogalaxy/imaging/fit_imaging.py
+++ b/autogalaxy/imaging/fit_imaging.py
@@ -7,7 +7,6 @@
 
 from autogalaxy.abstract_fit import AbstractFitInversion
 from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages
-from autogalaxy.analysis.preloads import Preloads
 from autogalaxy.galaxy.galaxy import Galaxy
 from autogalaxy.galaxy.galaxies import Galaxies
 from autogalaxy.galaxy.to_inversion import GalaxiesToInversion
@@ -26,7 +25,6 @@ def __init__(
         dataset_model: Optional[aa.DatasetModel] = None,
         adapt_images: Optional[AdaptImages] = None,
         settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads: aa.Preloads = Preloads(),
         run_time_dict: Optional[Dict] = None,
     ):
         """
@@ -66,16 +64,12 @@ def __init__(
             reconstructed galaxy's morphology.
         settings_inversion
             Settings controlling how an inversion is fitted for example which linear algebra formalism is used.
-        preloads
-            Contains preloaded calculations (e.g. linear algebra matrices) which can skip certain calculations in
-            the fit.
         run_time_dict
             A dictionary which if passed to the fit records how long fucntion calls which have the `profile_func`
             decorator take to run.
         """
 
         self.galaxies = Galaxies(galaxies=galaxies, run_time_dict=run_time_dict)
-        self.preloads = preloads
 
         super().__init__(
             dataset=dataset,
@@ -135,7 +129,6 @@ def galaxies_to_inversion(self) -> GalaxiesToInversion:
             galaxies=self.galaxies,
             adapt_images=self.adapt_images,
             settings_inversion=self.settings_inversion,
-            preloads=self.preloads,
             run_time_dict=self.run_time_dict,
         )
 
@@ -220,9 +213,9 @@ def subtracted_images_of_galaxies_dict(self) -> Dict[Galaxy, aa.Array2D]:
 
             subtracted_image = self.data - sum(other_galaxies_model_images)
 
-            subtracted_images_of_galaxies_dict[
-                self.galaxies[galaxy_index]
-            ] = subtracted_image
+            subtracted_images_of_galaxies_dict[self.galaxies[galaxy_index]] = (
+                subtracted_image
+            )
 
         return subtracted_images_of_galaxies_dict
 
@@ -285,44 +278,3 @@ def galaxies_linear_light_profiles_to_light_profiles(self) -> List[Galaxy]:
         or `GalaxyPlotter` objects.
         """
         return self.model_obj_linear_light_profiles_to_light_profiles
-
-    def refit_with_new_preloads(
-        self,
-        preloads: Preloads,
-        settings_inversion: Optional[aa.SettingsInversion] = None,
-    ) -> "FitImaging":
-        """
-        Returns a new fit which uses the dataset, galaxies and other objects of this fit, but uses a different set of
-        preloads input into this function.
-
-        This is used when setting up the preloads objects, to concisely test how using different preloads objects
-        changes the attributes of the fit.
-
-        Parameters
-        ----------
-        preloads
-            The new preloads which are used to refit the data using the
-        settings_inversion
-            Settings controlling how an inversion is fitted for example which linear algebra formalism is used.
-
-        Returns
-        -------
-        A new fit which has used new preloads input into this function but the same dataset, galaxies and other settings.
-        """
-        run_time_dict = {} if self.run_time_dict is not None else None
-
-        settings_inversion = (
-            self.settings_inversion
-            if settings_inversion is None
-            else settings_inversion
-        )
-
-        return FitImaging(
-            dataset=self.dataset,
-            galaxies=self.galaxies,
-            dataset_model=self.dataset_model,
-            adapt_images=self.adapt_images,
-            settings_inversion=settings_inversion,
-            preloads=preloads,
-            run_time_dict=run_time_dict,
-        )
diff --git a/autogalaxy/imaging/model/analysis.py b/autogalaxy/imaging/model/analysis.py
index 82f34282..ad2bc643 100644
--- a/autogalaxy/imaging/model/analysis.py
+++ b/autogalaxy/imaging/model/analysis.py
@@ -9,7 +9,6 @@
 
 from autogalaxy.analysis.adapt_images.adapt_image_maker import AdaptImageMaker
 from autogalaxy.analysis.analysis.dataset import AnalysisDataset
-from autogalaxy.analysis.preloads import Preloads
 from autogalaxy.cosmology.lensing import LensingCosmology
 from autogalaxy.cosmology.wrap import Planck15
 from autogalaxy.imaging.model.result import ResultImaging
@@ -94,9 +93,6 @@ def modify_before_fit(self, paths: af.DirectoryPaths, model: af.Collection):
 
         super().modify_before_fit(paths=paths, model=model)
 
-        if not paths.is_complete:
-            self.set_preloads(paths=paths, model=model)
-
         return self
 
     def log_likelihood_function(self, instance: af.ModelInstance) -> float:
@@ -153,7 +149,6 @@ def log_likelihood_function(self, instance: af.ModelInstance) -> float:
     def fit_from(
         self,
         instance: af.ModelInstance,
-        preload_overwrite: Optional[Preloads] = None,
         run_time_dict: Optional[Dict] = None,
     ) -> FitImaging:
         """
@@ -186,15 +181,12 @@ def fit_from(
 
         adapt_images = self.adapt_images_via_instance_from(instance=instance)
 
-        preloads = self.preloads if preload_overwrite is None else preload_overwrite
-
         return FitImaging(
             dataset=self.dataset,
             galaxies=galaxies,
             dataset_model=dataset_model,
             adapt_images=adapt_images,
             settings_inversion=self.settings_inversion,
-            preloads=preloads,
             run_time_dict=run_time_dict,
         )
 
@@ -253,7 +245,7 @@ def profile_log_likelihood_function(
         An `info_dict` is also created which stores information on aspects of the model and dataset that dictate
         run times, so the profiled times can be interpreted with this context.
 
-        The results of this profiling are then output to hard-disk in the `preloads` folder of the model-fit results,
+        The results of this profiling are then output to hard-disk in the `profiling` folder of the model-fit results,
         which they can be inspected to ensure run-times are as expected.
 
         Parameters
diff --git a/autogalaxy/interferometer/fit_interferometer.py b/autogalaxy/interferometer/fit_interferometer.py
index 7f7f6fe6..7095ced4 100644
--- a/autogalaxy/interferometer/fit_interferometer.py
+++ b/autogalaxy/interferometer/fit_interferometer.py
@@ -7,7 +7,6 @@
 
 from autogalaxy.abstract_fit import AbstractFitInversion
 from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages
-from autogalaxy.analysis.preloads import Preloads
 from autogalaxy.galaxy.galaxy import Galaxy
 from autogalaxy.galaxy.galaxies import Galaxies
 from autogalaxy.galaxy.to_inversion import GalaxiesToInversion
@@ -21,7 +20,6 @@ def __init__(
         dataset_model: Optional[aa.DatasetModel] = None,
         adapt_images: Optional[AdaptImages] = None,
         settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads: aa.Preloads = Preloads(),
         run_time_dict: Optional[Dict] = None,
     ):
         """
@@ -62,9 +60,6 @@ def __init__(
             reconstructed galaxy's morphology.
         settings_inversion
             Settings controlling how an inversion is fitted for example which linear algebra formalism is used.
-        preloads
-            Contains preloaded calculations (e.g. linear algebra matrices) which can skip certain calculations in
-            the fit.
         run_time_dict
             A dictionary which if passed to the fit records how long fucntion calls which have the `profile_func`
             decorator take to run.
@@ -92,8 +87,6 @@ def __init__(
         self.adapt_images = adapt_images
         self.settings_inversion = settings_inversion
 
-        self.preloads = preloads
-
     @property
     def profile_visibilities(self) -> aa.Visibilities:
         """
@@ -126,7 +119,6 @@ def galaxies_to_inversion(self) -> GalaxiesToInversion:
             galaxies=self.galaxies,
             adapt_images=self.adapt_images,
             settings_inversion=self.settings_inversion,
-            preloads=self.preloads,
             run_time_dict=self.run_time_dict,
         )
 
@@ -223,44 +215,3 @@ def galaxies_linear_light_profiles_to_light_profiles(self) -> List[Galaxy]:
         or `GalaxyPlotter` objects.
         """
         return self.model_obj_linear_light_profiles_to_light_profiles
-
-    def refit_with_new_preloads(
-        self,
-        preloads: Preloads,
-        settings_inversion: Optional[aa.SettingsInversion] = None,
-    ) -> "FitInterferometer":
-        """
-        Returns a new fit which uses the dataset, galaxies and other objects of this fit, but uses a different set of
-        preloads input into this function.
-
-        This is used when setting up the preloads objects, to concisely test how using different preloads objects
-        changes the attributes of the fit.
-
-        Parameters
-        ----------
-        preloads
-            The new preloads which are used to refit the data using the
-        settings_inversion
-            Settings controlling how an inversion is fitted for example which linear algebra formalism is used.
-
-        Returns
-        -------
-        A new fit which has used new preloads input into this function but the same dataset, galaxies and other settings.
-        """
-        if self.run_time_dict is not None:
-            run_time_dict = {}
-        else:
-            run_time_dict = None
-
-        if settings_inversion is None:
-            settings_inversion = self.settings_inversion
-
-        return FitInterferometer(
-            dataset=self.dataset,
-            galaxies=self.galaxies,
-            dataset_model=self.dataset_model,
-            adapt_images=self.adapt_images,
-            settings_inversion=settings_inversion,
-            preloads=preloads,
-            run_time_dict=run_time_dict,
-        )
diff --git a/autogalaxy/interferometer/model/analysis.py b/autogalaxy/interferometer/model/analysis.py
index c152e4e4..7ee73d1f 100644
--- a/autogalaxy/interferometer/model/analysis.py
+++ b/autogalaxy/interferometer/model/analysis.py
@@ -11,7 +11,6 @@
 
 from autogalaxy.analysis.adapt_images.adapt_image_maker import AdaptImageMaker
 from autogalaxy.analysis.analysis.dataset import AnalysisDataset
-from autogalaxy.analysis.preloads import Preloads
 from autogalaxy.cosmology.lensing import LensingCosmology
 from autogalaxy.cosmology.wrap import Planck15
 from autogalaxy.interferometer.model.result import ResultInterferometer
@@ -99,13 +98,6 @@ def modify_before_fit(self, paths: af.DirectoryPaths, model: af.Collection):
 
         super().modify_before_fit(paths=paths, model=model)
 
-        if not paths.is_complete:
-            logger.info(
-                "PRELOADS - Setting up preloads, may take a few minutes for fits using an inversion."
-            )
-
-            self.set_preloads(paths=paths, model=model)
-
         return self
 
     def log_likelihood_function(self, instance: af.ModelInstance) -> float:
@@ -161,7 +153,6 @@ def log_likelihood_function(self, instance: af.ModelInstance) -> float:
     def fit_from(
         self,
         instance: af.ModelInstance,
-        preload_overwrite: Optional[Preloads] = None,
         run_time_dict: Optional[Dict] = None,
     ) -> FitInterferometer:
         """
@@ -191,14 +182,11 @@ def fit_from(
 
         adapt_images = self.adapt_images_via_instance_from(instance=instance)
 
-        preloads = self.preloads if preload_overwrite is None else preload_overwrite
-
         return FitInterferometer(
             dataset=self.dataset,
             galaxies=galaxies,
             adapt_images=adapt_images,
             settings_inversion=self.settings_inversion,
-            preloads=preloads,
             run_time_dict=run_time_dict,
         )
 
@@ -259,7 +247,7 @@ def profile_log_likelihood_function(
         An `info_dict` is also created which stores information on aspects of the model and dataset that dictate
         run times, so the profiled times can be interpreted with this context.
 
-        The results of this profiling are then output to hard-disk in the `preloads` folder of the model-fit results,
+        The results of this profiling are then output to hard-disk in the `profiling` folder of the model-fit results,
         which they can be inspected to ensure run-times are as expected.
 
         Parameters
diff --git a/autogalaxy/operate/deflections.py b/autogalaxy/operate/deflections.py
index 28f01958..57af2276 100644
--- a/autogalaxy/operate/deflections.py
+++ b/autogalaxy/operate/deflections.py
@@ -662,9 +662,7 @@ def einstein_mass_angular_list_from(
         einstein_radius_list = self.einstein_radius_list_from(
             grid=grid, pixel_scale=pixel_scale
         )
-        return [
-            np.pi * einstein_radius**2 for einstein_radius in einstein_radius_list
-        ]
+        return [np.pi * einstein_radius**2 for einstein_radius in einstein_radius_list]
 
     @evaluation_grid
     def einstein_mass_angular_from(
diff --git a/autogalaxy/profiles/basis.py b/autogalaxy/profiles/basis.py
index 63747b22..f71777a0 100644
--- a/autogalaxy/profiles/basis.py
+++ b/autogalaxy/profiles/basis.py
@@ -125,9 +125,11 @@ def image_2d_list_from(
         The image of the light profiles in the basis summed together.
         """
         return [
-            light_profile.image_2d_from(grid=grid, operated_only=operated_only)
-            if not isinstance(light_profile, lp_linear.LightProfileLinear)
-            else np.zeros((grid.shape[0],))
+            (
+                light_profile.image_2d_from(grid=grid, operated_only=operated_only)
+                if not isinstance(light_profile, lp_linear.LightProfileLinear)
+                else np.zeros((grid.shape[0],))
+            )
             for light_profile in self.light_profile_list
         ]
 
diff --git a/autogalaxy/profiles/geometry_profiles.py b/autogalaxy/profiles/geometry_profiles.py
index 725a1049..64fb0c45 100644
--- a/autogalaxy/profiles/geometry_profiles.py
+++ b/autogalaxy/profiles/geometry_profiles.py
@@ -311,7 +311,8 @@ def elliptical_radii_grid_from(
         """
         return np.sqrt(
             np.add(
-                np.square(grid.array[:, 1]), np.square(np.divide(grid.array[:, 0], self.axis_ratio))
+                np.square(grid.array[:, 1]),
+                np.square(np.divide(grid.array[:, 0], self.axis_ratio)),
             )
         )
 
@@ -336,7 +337,7 @@ def eccentric_radii_grid_from(
 
         grid_radii = self.elliptical_radii_grid_from(grid=grid, **kwargs).array
 
-        return np.multiply(np.sqrt(self.axis_ratio), grid_radii)#.view(np.ndarray)
+        return np.multiply(np.sqrt(self.axis_ratio), grid_radii)  # .view(np.ndarray)
 
     @aa.grid_dec.to_grid
     def transformed_to_reference_frame_grid_from(
diff --git a/autogalaxy/profiles/mass/abstract/jax_utils.py b/autogalaxy/profiles/mass/abstract/jax_utils.py
index 735342ef..dd31478c 100644
--- a/autogalaxy/profiles/mass/abstract/jax_utils.py
+++ b/autogalaxy/profiles/mass/abstract/jax_utils.py
@@ -26,7 +26,7 @@ def reg4(z, _, i_sqrt_pi):
 
 
 def reg5(z, sqrt_pi, _):
-    mz2 = -z**2
+    mz2 = -(z**2)
     f1 = sqrt_pi
     f2 = 1.0
     s1 = [1.320522, 35.7668, 219.031, 1540.787, 3321.99, 36183.31]
@@ -80,7 +80,7 @@ def w_f_approx(z):
     r2 = (abs_z2 >= 256.0) & (abs_z2 < 38000.0)
     r3 = (abs_z2 >= 62.0) & (abs_z2 < 256.0)
     r4 = (abs_z2 >= 30.0) & (abs_z2 < 62.0) & (z_imag2 >= 1e-13)
-    # region bounds for 5 taken directly from Zaghloul (2017) 
+    # region bounds for 5 taken directly from Zaghloul (2017)
     # https://dl.acm.org/doi/pdf/10.1145/3119904
     r5_1 = (abs_z2 >= 30.0) & (abs_z2 < 62.0) & (z_imag2 < 1e-13)
     r5_2 = (abs_z2 >= 2.5) & (abs_z2 < 30.0) & (z_imag2 < 0.072)
@@ -103,8 +103,8 @@ def w_f_approx_jvp(primals, tangents):
     # define a custom jvp to avoid the issue using `jnp.where` with `jax.grad`
     # also the derivative is defined analytically for this function so bypass
     # auto diffing over the complex functions above.
-    z, = primals
-    z_dot, = tangents
+    (z,) = primals
+    (z_dot,) = tangents
     primal_out = w_f_approx(z)
     i_sqrt_pi = 1j / jnp.sqrt(jnp.pi)
     tangent_out = z_dot * 2 * (i_sqrt_pi - z * primal_out)
diff --git a/autogalaxy/profiles/mass/abstract/mge_jax.py b/autogalaxy/profiles/mass/abstract/mge_jax.py
index a9845258..2fbfd5ed 100644
--- a/autogalaxy/profiles/mass/abstract/mge_jax.py
+++ b/autogalaxy/profiles/mass/abstract/mge_jax.py
@@ -65,16 +65,13 @@ def eta(p):
         """
 
         i = np.arange(1, p, 1)
-        kesi_last = 1/2**p
+        kesi_last = 1 / 2**p
         k = kesi_last + np.cumsum(np.cumprod((p + 1 - i) / i) * kesi_last)
-        
-        kesi_list = np.hstack([
-            np.array([0.5]),
-            np.ones(p),
-            k[::-1],
-            np.array([kesi_last])
-        ])
-        coef = (-1)**np.arange(0, 2 * p + 1, 1)
+
+        kesi_list = np.hstack(
+            [np.array([0.5]), np.ones(p), k[::-1], np.array([kesi_last])]
+        )
+        coef = (-1) ** np.arange(0, 2 * p + 1, 1)
         eta_list = coef * 2.0 * np.sqrt(2.0 * np.pi) * 10 ** (p / 3.0) * kesi_list
         return eta_list
 
@@ -111,7 +108,9 @@ def _decompose_convergence_via_mge(
         sigma_list = np.exp(log_sigmas)
 
         amplitude_list = np.zeros(func_gaussians)
-        f_sigma = np.sum(etas * np.real(func(sigma_list.reshape(-1, 1) * kesis)), axis=1)
+        f_sigma = np.sum(
+            etas * np.real(func(sigma_list.reshape(-1, 1) * kesis)), axis=1
+        )
         amplitude_list = f_sigma * d_log_sigma / np.sqrt(2.0 * np.pi)
         amplitude_list = amplitude_list.at[0].multiply(0.5)
         amplitude_list = amplitude_list.at[-1].multiply(0.5)
@@ -120,7 +119,9 @@ def _decompose_convergence_via_mge(
     def convergence_2d_via_mge_from(self, grid_radii):
         raise NotImplementedError()
 
-    def _convergence_2d_via_mge_from(self, grid_radii, func_terms=28, func_gaussians=20):
+    def _convergence_2d_via_mge_from(
+        self, grid_radii, func_terms=28, func_gaussians=20
+    ):
         """Calculate the projected convergence at a given set of arc-second gridded coordinates.
 
         Parameters
@@ -129,19 +130,25 @@ def _convergence_2d_via_mge_from(self, grid_radii, func_terms=28, func_gaussians
             The grid of (y,x) arc-second coordinates the convergence is computed on.
 
         """
-        amps, sigmas = self.decompose_convergence_via_mge(func_terms=func_terms, func_gaussians=func_gaussians)
+        amps, sigmas = self.decompose_convergence_via_mge(
+            func_terms=func_terms, func_gaussians=func_gaussians
+        )
 
         convergence = 0.0
 
         inv_sigma_ = 1 / sigmas.reshape((-1,) + (1,) * grid_radii.array.ndim)
         amps_ = amps.reshape((-1,) + (1,) * grid_radii.array.ndim)
-        convergence = amps_ * np.exp(-0.5 * (grid_radii.array * inv_sigma_)**2)
+        convergence = amps_ * np.exp(-0.5 * (grid_radii.array * inv_sigma_) ** 2)
         return convergence.sum(axis=0)
 
-    def _deflections_2d_via_mge_from(self, grid, sigmas_factor=1.0, func_terms=28, func_gaussians=20):
+    def _deflections_2d_via_mge_from(
+        self, grid, sigmas_factor=1.0, func_terms=28, func_gaussians=20
+    ):
         axis_ratio = np.min(np.array([self.axis_ratio, 0.9999]))
 
-        amps, sigmas = self.decompose_convergence_via_mge(func_terms=func_terms, func_gaussians=func_gaussians)
+        amps, sigmas = self.decompose_convergence_via_mge(
+            func_terms=func_terms, func_gaussians=func_gaussians
+        )
         sigmas *= sigmas_factor
 
         angle = self.zeta_from(
diff --git a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
index 9247c9c9..838b6745 100644
--- a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
+++ b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
@@ -43,9 +43,7 @@ def kappa_s_and_scale_radius(
 
     ##############################
     def integrand(r):
-        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (
-            inner_slope - 3
-        )
+        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (inner_slope - 3)
 
     de_c = (
         (overdens / 3.0)
diff --git a/autogalaxy/profiles/mass/dark/nfw_hk24_util.py b/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
index 0c46b4cf..71bdce62 100644
--- a/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
+++ b/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
@@ -3,6 +3,7 @@
 
 @author: felixvecchi
 """
+
 import numpy as np
 from astropy.cosmology import Planck15
 
@@ -206,16 +207,11 @@ def g1_g2_from(x1, x2, e, k_s):
 
     # Prefactors for g1
     pre_f0_g1 = (
-        ((x1 - e) ** 2 + x2**2)
-        * ((x1 + e) ** 2 + x2**2)
-        * (x1**2 - x2**2 - e**2)
+        ((x1 - e) ** 2 + x2**2) * ((x1 + e) ** 2 + x2**2) * (x1**2 - x2**2 - e**2)
     )
 
     pre_f1_g1 = (
-        2
-        * e**2
-        * (x1**2 - 1)
-        * ((x1**2 - x2**2 - e**2) ** 2 - 4 * x1**2 * x2**2)
+        2 * e**2 * (x1**2 - 1) * ((x1**2 - x2**2 - e**2) ** 2 - 4 * x1**2 * x2**2)
         - (x1**2 - x2**2 - e**2) ** 3 * (3 + e**2) / 2
         + 6 * x1**2 * x2**2 * (e**2 - 1) * (x1**2 - x2**2 - e**2)
     )
@@ -225,12 +221,7 @@ def g1_g2_from(x1, x2, e, k_s):
         - 8 * e**2 * x1**2 * x2**2
     )
 
-    pre_f3_g1 = (
-        2
-        * x1
-        * x2
-        * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
-    )
+    pre_f3_g1 = 2 * x1 * x2 * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
 
     # First component shear
     g1 = full_pre_factor * (
@@ -245,21 +236,12 @@ def g1_g2_from(x1, x2, e, k_s):
         * x2
         * (
             (x1**2 + x2**2 + e**2)
-            * (
-                (5 * e**2 - 3) * x1**2
-                - 3 * (1 + e**2) * x2**2
-                + (5 - 3 * e**2) * e**2
-            )
+            * ((5 * e**2 - 3) * x1**2 - 3 * (1 + e**2) * x2**2 + (5 - 3 * e**2) * e**2)
             - 4 * e**2 * x1**2 * (1 + e**2)
         )
     )
 
-    pre_f2_g2 = (
-        -2
-        * x1
-        * x2
-        * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
-    )
+    pre_f2_g2 = -2 * x1 * x2 * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
 
     pre_f3_g2 = -(
         (x1**2 - x2**2 - e**2) * ((x1**2 + x2**2) ** 2 - e**4)
diff --git a/autogalaxy/profiles/mass/stellar/sersic.py b/autogalaxy/profiles/mass/stellar/sersic.py
index 548a5166..a7608c90 100644
--- a/autogalaxy/profiles/mass/stellar/sersic.py
+++ b/autogalaxy/profiles/mass/stellar/sersic.py
@@ -138,7 +138,9 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs):
     @aa.grid_dec.to_vector_yx
     @aa.grid_dec.transform
     @aa.grid_dec.relocate_to_radial_minimum
-    def deflections_2d_via_mge_from(self, grid: aa.type.Grid2DLike, func_terms=28, func_gaussians=20, **kwargs):
+    def deflections_2d_via_mge_from(
+        self, grid: aa.type.Grid2DLike, func_terms=28, func_gaussians=20, **kwargs
+    ):
         """
         Calculate the projected 2D deflection angles from a grid of (y,x) arc second coordinates, by computing and
         summing the convergence of each individual cse used to decompose the mass profile.
@@ -153,7 +155,10 @@ def deflections_2d_via_mge_from(self, grid: aa.type.Grid2DLike, func_terms=28, f
             The grid of (y,x) arc-second coordinates the convergence is computed on.
         """
         return self._deflections_2d_via_mge_from(
-            grid=grid, sigmas_factor=np.sqrt(self.axis_ratio), func_terms=func_terms, func_gaussians=func_gaussians,
+            grid=grid,
+            sigmas_factor=np.sqrt(self.axis_ratio),
+            func_terms=func_terms,
+            func_gaussians=func_gaussians,
         )
 
     @aa.grid_dec.to_vector_yx
@@ -196,7 +201,9 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, **kwargs):
     @aa.grid_dec.to_array
     @aa.grid_dec.transform
     @aa.grid_dec.relocate_to_radial_minimum
-    def convergence_2d_via_mge_from(self, grid: aa.type.Grid2DLike, func_terms=28, func_gaussians=20, **kwargs):
+    def convergence_2d_via_mge_from(
+        self, grid: aa.type.Grid2DLike, func_terms=28, func_gaussians=20, **kwargs
+    ):
         """
         Calculate the projected convergence at a given set of arc-second gridded coordinates.
 
@@ -209,7 +216,11 @@ def convergence_2d_via_mge_from(self, grid: aa.type.Grid2DLike, func_terms=28, f
 
         eccentric_radii = self.eccentric_radii_grid_from(grid=grid, **kwargs)
 
-        return self._convergence_2d_via_mge_from(grid_radii=eccentric_radii, func_terms=func_terms, func_gaussians=func_gaussians)
+        return self._convergence_2d_via_mge_from(
+            grid_radii=eccentric_radii,
+            func_terms=func_terms,
+            func_gaussians=func_gaussians,
+        )
 
     @aa.over_sample
     @aa.grid_dec.to_array
@@ -255,7 +266,9 @@ def image_2d_via_radii_from(self, radius: np.ndarray):
             * (((radius / self.effective_radius) ** (1.0 / self.sersic_index)) - 1)
         )
 
-    def decompose_convergence_via_mge(self, func_terms=28, func_gaussians=20) -> Tuple[List, List]:
+    def decompose_convergence_via_mge(
+        self, func_terms=28, func_gaussians=20
+    ) -> Tuple[List, List]:
         radii_min = self.effective_radius / 100.0
         radii_max = self.effective_radius * 20.0
 
@@ -270,7 +283,11 @@ def sersic_2d(r):
             )
 
         return self._decompose_convergence_via_mge(
-            func=sersic_2d, radii_min=radii_min, radii_max=radii_max, func_terms=func_terms, func_gaussians=func_gaussians
+            func=sersic_2d,
+            radii_min=radii_min,
+            radii_max=radii_max,
+            func_terms=func_terms,
+            func_gaussians=func_gaussians,
         )
 
     def decompose_convergence_via_cse(
diff --git a/autogalaxy/profiles/mass/total/isothermal.py b/autogalaxy/profiles/mass/total/isothermal.py
index 734330eb..79f37893 100644
--- a/autogalaxy/profiles/mass/total/isothermal.py
+++ b/autogalaxy/profiles/mass/total/isothermal.py
@@ -40,10 +40,12 @@ def psi_from(grid, axis_ratio, core_radius):
     """
     if USING_JAX:
         return np.sqrt(
-            (axis_ratio**2.0 * (grid[:, 1]**2.0 + core_radius**2.0)) + grid[:, 0]**2.0 + 1e-16
+            (axis_ratio**2.0 * (grid[:, 1] ** 2.0 + core_radius**2.0))
+            + grid[:, 0] ** 2.0
+            + 1e-16
         )
     else:
-     return np.sqrt(
+        return np.sqrt(
             np.add(
                 np.multiply(
                     axis_ratio**2.0, np.add(np.square(grid[:, 1]), core_radius**2.0)
diff --git a/autogalaxy/profiles/mass/total/jax_utils.py b/autogalaxy/profiles/mass/total/jax_utils.py
index 42fa0026..5bed36c6 100644
--- a/autogalaxy/profiles/mass/total/jax_utils.py
+++ b/autogalaxy/profiles/mass/total/jax_utils.py
@@ -5,7 +5,7 @@
 
 # A version of scan that will *not* re-compile partial functions when variables change
 # taken from https://github.com/google/jax/issues/14743#issuecomment-1456900634
-scan = jax.jit(jax.lax.scan, static_argnames=('length', 'reverse', 'unroll'))
+scan = jax.jit(jax.lax.scan, static_argnames=("length", "reverse", "unroll"))
 
 
 def body_fun(carry, n, factor, ei2phi, slope):
@@ -19,7 +19,7 @@ def body_fun(carry, n, factor, ei2phi, slope):
 
 
 def omega(eiphi, slope, factor, n_terms=20):
-    '''JAX implementation of the numerical evaluation of the angular component of
+    """JAX implementation of the numerical evaluation of the angular component of
     the complex deflection angle for the elliptical power law profile as given as
     given by Tessore and Metcalf 2015.  Based on equation 29, and gives
     omega (e.g. can be used as a drop in replacement for the exp(i * phi) * special.hyp2f1
@@ -37,11 +37,11 @@ def omega(eiphi, slope, factor, n_terms=20):
     n_terms:
         The number of terms to calculate for the series expansion, defaults to 20 (this should
         be sufficient most of the time)
-    '''
+    """
     # use modified scan with a partial'ed function to avoid re-compile
     (_, partial_sum), _ = scan(
         partial(body_fun, factor=factor, ei2phi=eiphi**2, slope=slope),
         (eiphi, eiphi),
-        jnp.arange(1, n_terms)
+        jnp.arange(1, n_terms),
     )
     return partial_sum
diff --git a/autogalaxy/profiles/mass/total/power_law.py b/autogalaxy/profiles/mass/total/power_law.py
index 99404cfa..cdab5afe 100644
--- a/autogalaxy/profiles/mass/total/power_law.py
+++ b/autogalaxy/profiles/mass/total/power_law.py
@@ -91,25 +91,21 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, **kwargs):
         z = np.add(
             np.multiply(np.cos(angle), 1 + 0j), np.multiply(np.sin(angle), 0 + 1j)
         )
-       
+
         if USING_JAX:
             # offset radius so calculation is finite at (0, 0)
-            R = np.sqrt(
-                (self.axis_ratio * grid[:, 1])**2 + grid[:, 0]**2 + 1e-16
-            )
+            R = np.sqrt((self.axis_ratio * grid[:, 1]) ** 2 + grid[:, 0] ** 2 + 1e-16)
             zh = omega(z, slope, factor, n_terms=20)
         else:
             R = np.sqrt(
-                np.add(np.multiply(self.axis_ratio**2, grid[:, 1] ** 2), grid[:, 0] ** 2)
+                np.add(
+                    np.multiply(self.axis_ratio**2, grid[:, 1] ** 2), grid[:, 0] ** 2
+                )
             )
             zh = z * special.hyp2f1(1.0, 0.5 * slope, 2.0 - 0.5 * slope, -factor * z**2)
 
         complex_angle = (
-            2.0
-            * b
-            / (1.0 + self.axis_ratio)
-            * (b / R) ** (slope - 1.0)
-            * zh
+            2.0 * b / (1.0 + self.axis_ratio) * (b / R) ** (slope - 1.0) * zh
         )
 
         deflection_y = complex_angle.imag
diff --git a/autogalaxy/util/shear_field.py b/autogalaxy/util/shear_field.py
index b0f5f02d..179e8f09 100644
--- a/autogalaxy/util/shear_field.py
+++ b/autogalaxy/util/shear_field.py
@@ -94,7 +94,6 @@ class ShearYX2D(aa.VectorYX2D, AbstractShearField):
 
 
 class ShearYX2DIrregular(aa.VectorYX2DIrregular, AbstractShearField):
-
     """
     An irregular shear field, which is collection of (y,x) vectors which are located on an irregular grid
     of (y,x) coordinates.
diff --git a/docs/api/plot.rst b/docs/api/plot.rst
index 33219479..3fe971ef 100644
--- a/docs/api/plot.rst
+++ b/docs/api/plot.rst
@@ -83,7 +83,7 @@ visuals to figures.
 Matplot Lib Wrappers [aplt]
 ---------------------------
 
-Wrappers for every ``matplotlib`` function used by a ``Plotter``, allowing for detailed customizaiton of
+Wrappers for every ``matplotlib`` function used by a ``Plotter``, allowing for detailed customization of
 every figure and subplot.
 
 .. currentmodule:: autogalaxy.plot
diff --git a/docs/installation/conda.rst b/docs/installation/conda.rst
index 27f222bf..6957e662 100644
--- a/docs/installation/conda.rst
+++ b/docs/installation/conda.rst
@@ -46,7 +46,7 @@ You may get warnings which state something like:
 
 .. code-block:: bash
 
-    ERROR: autoarray 2024.5.16.0 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible.
+    ERROR: autoarray 2024.9.21.2 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible.
     ERROR: numba 0.53.1 has requirement llvmlite<0.37,>=0.36.0rc1, but you'll have llvmlite 0.38.0 which is incompatible.
 
 If you see these messages, they do not mean that the installation has failed and the instructions below will
diff --git a/docs/installation/pip.rst b/docs/installation/pip.rst
index c8321746..760f8667 100644
--- a/docs/installation/pip.rst
+++ b/docs/installation/pip.rst
@@ -27,7 +27,7 @@ You may get warnings which state something like:
 
 .. code-block:: bash
 
-    ERROR: autoarray 2024.5.16.0 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible.
+    ERROR: autoarray 2024.9.21.2 has requirement numpy<=1.22.1, but you'll have numpy 1.22.2 which is incompatible.
     ERROR: numba 0.53.1 has requirement llvmlite<0.37,>=0.36.0rc1, but you'll have llvmlite 0.38.0 which is incompatible.
 
 If you see these messages, they do not mean that the installation has failed and the instructions below will
diff --git a/docs/overview/overview_2_new_user_guide.rst b/docs/overview/overview_2_new_user_guide.rst
index 91614f82..8904e200 100644
--- a/docs/overview/overview_2_new_user_guide.rst
+++ b/docs/overview/overview_2_new_user_guide.rst
@@ -3,28 +3,241 @@
 New User Guide
 ==============
 
-1) Read the `autolens_workspace/start_here.ipynb` notebook if yo have not already.
+**PyAutoGalaxy** is an extensive piece of software with functionality for doing many different analysis tasks, fitting
+different data types and it is used for a variety of different science cases. This means the documentation is quite
+extensive, and it may be difficult to find the example script you need.
 
-2) Decide if you should take the HowToLens lectures.
+This page provides a sequential guide for news users on how to begin learning **PyAutoGalaxy**, and can act as a useful
+resource for existing users who are looking for how to do a specific task.
 
-3) Config files and make sure default visualization settings are correct.
+Before starting this guide, you should ensure you have installed **PyAutoGalaxy** and downloaded the ``autogalaxy_workspace``
+by following the `installation guide <https://PyAutoGalaxy.readthedocs.io/en/latest/installation/overview.html>`_.
 
-4) Determine if you are `imaging` / `interferometer` / `point_source` / `group`.
+Contents
+--------
 
-5) Quickly look at the API reference guides and unit convention guides, so you know to refer to them if you get stuck (`autolens_workspace/api`).
+One line summaries of each step in the new user guide is given below, to give you a sense of what you are going to learn:
 
-6) Learn how to simulate data (`autolens_workspace/start_here/simulators/start_here.ipynb`).
+**1) Workspace:** Read the ``start_here.ipynb`` workspace example for a quick run through of the core API.
+**2) HowToGalaxy?**: Whether you should begin with lectures aimed at inexperienced scientists (e.g. under graduate students).
 
-7) Learn how to model data (`autolens_workspace/start_here/modeling/start_here.ipynb`).
 
-8) Learn how to prepare data and model it (`autolens_workspace/start_here/data_preparation/start_here.ipynb`).
+1) Workspace
+------------
 
-9) Learn more about results (`autolens_workspace/start_here/results/start_here.ipynb`).
+You should now have the ``autogalaxy_workspace`` on your computer and see many of the folder and files we'll begin
+navigating.
 
-10) Learn how to plot data (`autolens_workspace/start_here/plotting/start_here.ipynb`).
+First of all, if you have not already, you should read the `autogalaxy_workspace/start_here.ipynb` notebook,
+which provides a run through of the core API for galaxy morphology calculations and modeling.
 
-11) Decide if you need to use any number of features in the `autolens_workspace/features` package.
+GitHub Links:
 
-12) Advanced stuff.
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release
 
-13) Do science and cite correctly!
\ No newline at end of file
+2) HowToGalaxy?
+-------------
+
+For experienced scientists, the **PyAutoGalaxy** examples will be simple to follow. Concepts surrounding galaxy morphology may
+already be familiar and the statistical techniques used for fitting and modeling already understood.
+
+For those less familiar with these concepts (e.g. undergraduate students, new PhD students or interested members of the
+public), things may have been less clear and a slower more detailed explanation of each concept would be beneficial.
+
+The **HowToGalaxy** Jupyter Notebook lectures provide exactly this. They are a 3+ chapter guide which thoroughly
+take you through the core concepts of galaxy morphology, teach you the principles of the statistical techniques
+used in modeling and ultimately will allow you to undertake scientific research like a professional astronomer.
+
+To complete thoroughly, they'll probably take 2-4 days, so you may want try moving ahead to the examples but can
+go back to these lectures if you find them hard to follow.
+
+If this sounds like it suits you, checkout the ``autogalaxy_workspace/notebooks/HowToGalaxy`` package now.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/HowToGalaxy
+
+3) Configs
+----------
+
+The ``autogalaxy_workspace/config`` folder contains numerous .YAML configuration files which customization many
+default settings of **PyAutoGalaxy**.
+
+Documentation for all config settings are provided within each config file.
+
+New users should not worry about the majority of configs for now. However, the ``config/visualize`` folder contains
+config files which customization ``matplotlib`` visualization, and editing these now will ensure figures and
+images display optimally in your Jupyter Notebooks.
+
+All default ``matplotlib`` options are customized via the `mat_wrap.yaml`, `mat_wrap_1d.yaml` and `mat_wrap_2d.yaml` files
+in `autogalaxy_workspace/config/visualize/mat_wrap`. For example, if figures display with labels that are too big
+or small, you can adjust their default labelsizes by changing the following options:
+
+ - mat_wrap.yaml -> Figure -> figure: -> figsize
+ - mat_wrap.yaml -> YLabel -> figure: -> fontsize
+ - mat_wrap.yaml -> XLabel -> figure: -> fontsize
+ - mat_wrap.yaml -> TickParams -> figure: -> labelsize
+ - mat_wrap.yaml -> YTicks -> figure: -> labelsize
+ - mat_wrap.yaml -> XTicks -> figure: -> labelsize
+
+The default colormap can be changed from the default to your favour ``matplotlib`` colormap, but adjusting:
+
+ - mat_wrap.yaml -> Cmap -> figure -> cmap
+
+All settings have a ``figure`` and ``subplot`` option, so that single image ``figures`` and a subplot of multiple
+figures can be customized independently.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/config
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/config/visualize
+
+4) Dataset Type
+---------------
+
+**PyAutoGalaxy** supports multiple different data types, and you as a user likely only require to learn how to use
+the software to analyse one type of dataset.
+
+Therefore, you now need to assess which dataset type is relevant to you:
+
+- **Imaging**: CCD imaging data (e.g. from the Hubble Space Telescope or James Webb Space Telescope), in which case
+you will go to the ``imaging`` packages in the workspace.
+
+- **Interferometry**: Interferometer data from a submm or radio interferometer (e.g. ALMA or JVLA), in which case
+you will go to the ``interferometer`` packages in the workspace.
+
+5) API and Units Guides
+-----------------------
+
+The ``autogalaxy_workspace/guides`` package has many useful guides, including concise API reference guides (``guides/api``)
+and unit conversion guides (``guides/units``).
+
+Quickly navigate to this part of the workspace and skim read the guides quickly. You do not need to understand them in detail now
+so don't spend long reading them.
+
+The purpose of looking at them now is you know they exist and can refer to them if you get stuck using **PyAutoGalaxy**.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/guides
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/guides/api
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/guides/units
+
+6) Simulations
+--------------
+
+Learning how to simulate your type of data is the best way to understanding how to analyse it.
+
+Therefore, in the ``autogalaxy_workspace/simulators`` folder, find the ``start_here.ipynb`` of your dataset.
+
+For example, if your dataset type is CCD imaging data, you'll read the notebook ``autogalaxy_workspace/simulators/imaging/start_here.ipynb``.
+
+Your **PyAutoGalaxy** use case might only require you to be able to simulate galaxies, for example if you are
+training a neural network. In this case, you can stop the guide and use the tools in the ``simulators`` package
+to start doing your science!
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/simulators
+
+7) Modeling
+-----------
+
+Having simulated a dataset, you are now ready to learn how to model it.
+
+Therefore, in the ``autogalaxy_workspace/modeling`` folder, find the ``start_here.ipynb`` of your dataset.
+
+For example, if your dataset type is CCD imaging data, you'll read the notebook ``autogalaxy_workspace/modeling/imaging/start_here.ipynb``.
+
+Your **PyAutoGalaxy** use case might only require you to be able to model simulated galaxies, for example if you are
+investigating what models can be used to learn about galaxy structure. In this case, you can skip the data preparation
+step below and go straight to learning about results.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/modeling
+
+8) Data Preparation
+-------------------
+
+If you have real observations of galaxies you want to model, you need to prepare the data so that it
+is appropriate for **PyAutoGalaxy**.
+
+This includes reducing the data so the galaxy is in the centre of the image, making sure all units
+are defined correctly and reducing extra data products like the Point Spread Function for CCD imaging data.
+
+Therefore, in the ``autogalaxy_workspace/data_preparation`` folder, find the ``start_here.ipynb`` of your dataset.
+
+For example, if your dataset type is CCD imaging data, you'll read the notebook ``autogalaxy_workspace/data_preparation/imaging/start_here.ipynb``.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/data_preparation
+
+9) Results
+----------
+
+Modeling infers many results, including parameter estimates, posteriors and a Bayesian evidence of the model.
+Furthermore, you may wish to inspect the results, the quality of the fit and produce visuals to determine
+if you think its a good fit.
+
+Therefore, now read the ``autogalaxy_workspace/*/results/start_here.ipynb`` notebook.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/results
+
+10) Plotting
+------------
+
+**PyAutoGalaxy** has an in depth visualizaiton library that allows for high levels of customization via ``matplotlib``.
+
+Plotting has its own dedicated API, which you should become familiar with via the example ``autogalaxy_workspace/*/plot/start_here.ipynb``.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/blob/main/notebooks/plot/start_here.ipynb
+
+11) Features
+------------
+
+You now have a comprehensive understanding of the **PyAutoGalaxy** API and how to use it to simulate, model and
+plot your data.
+
+**PyAutoGalaxy** has many more features, which may or may not be useful for your science case.
+
+Example notebooks for every feature are provided in the ``autogalaxy_workspace/*/features`` package and a high-level
+summary of each feature is provided on the next page of this readthedocs.
+
+What features you need depend on many factors: (i) your science case; (ii) the quality of your data; (iii) how
+much time you are willing to invest in learning **PyAutoGalaxy**. We recommend you read the literature in conjunction
+with assessing what features are available, and then make an informed decision on what is appropriate for you.
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/features
+
+12) Advanced
+------------
+
+The ``autogalaxy_workspace/*/advanced`` folder has numerous advanced examples which only a user experienced with
+**PyAutoGalaxy** should use.
+
+These include examples of how to fit multiple datasets simultaneously (e.g. multi-wavelength CCD imaging datasets),
+automated pipelines for modeling large galaxy samples (called the Source, Light and Mass (SLaM) pipelines in the
+literature) and a step-by-step guide of the **PyAutoGalaxy** likelihood function.
+
+New users should ignore this folder for now, but note that you may find it has important functionality for
+your science research in a couple of months time once you are experienced with **PyAutoGalaxy**!
+
+GitHub Links:
+
+https://github.com/Jammy2211/autogalaxy_workspace/tree/release/notebooks/advanced
+
+Wrap Up
+-------
+
+After completing this guide, you should be able to use **PyAutoGalaxy** for your science research.
+
+The biggest decisions you'll need to make are what features and functionality your specific science case requires,
+which the next readthedocs page gives an overview of to help you decide.
\ No newline at end of file


Test patch:
diff --git a/test_autogalaxy/analysis/test_preloads.py b/test_autogalaxy/analysis/test_preloads.py
deleted file mode 100644
index 28496242..00000000
--- a/test_autogalaxy/analysis/test_preloads.py
+++ /dev/null
@@ -1,130 +0,0 @@
-import numpy as np
-from os import path
-
-import autofit as af
-
-import autogalaxy as ag
-
-
-def test__set_blurred_image():
-    # Blurred image is all zeros so preloads as zeros
-
-    fit_0 = ag.m.MockFitImaging(blurred_image=np.zeros(2))
-    fit_1 = ag.m.MockFitImaging(blurred_image=np.zeros(2))
-
-    preloads = ag.Preloads(blurred_image=1)
-    preloads.set_blurred_image(fit_0=fit_0, fit_1=fit_1)
-
-    assert preloads.blurred_image is None
-
-    # Blurred image are different, indicating the model parameters change the grid, so no preloading.
-
-    fit_0 = ag.m.MockFitImaging(blurred_image=np.array([1.0]))
-    fit_1 = ag.m.MockFitImaging(blurred_image=np.array([2.0]))
-
-    preloads = ag.Preloads(blurred_image=1)
-    preloads.set_blurred_image(fit_0=fit_0, fit_1=fit_1)
-
-    assert preloads.blurred_image is None
-
-    # Blurred images are the same meaning they are fixed in the model, so do preload.
-
-    fit_0 = ag.m.MockFitImaging(blurred_image=np.array([1.0]))
-    fit_1 = ag.m.MockFitImaging(blurred_image=np.array([1.0]))
-
-    preloads = ag.Preloads(blurred_image=1)
-    preloads.set_blurred_image(fit_0=fit_0, fit_1=fit_1)
-
-    assert (preloads.blurred_image == np.array([1.0])).all()
-
-
-def test__info():
-    file_path = path.join("{}".format(path.dirname(path.realpath(__file__))), "files")
-
-    file_preloads = path.join(file_path, "preloads.summary")
-
-    preloads = ag.Preloads(
-        blurred_image=np.zeros(3),
-        w_tilde=None,
-        use_w_tilde=False,
-        image_plane_mesh_grid_pg_list=None,
-        relocated_grid=None,
-        mapper_list=None,
-        operated_mapping_matrix=None,
-    )
-
-    af.formatter.output_list_of_strings_to_file(
-        file=file_preloads, list_of_strings=preloads.info
-    )
-
-    results = open(file_preloads)
-    lines = results.readlines()
-
-    i = 0
-
-    assert lines[i] == f"W Tilde = False\n"
-    i += 1
-    assert lines[i] == f"Use W Tilde = False\n"
-    i += 1
-    assert lines[i] == f"\n"
-    i += 1
-    assert lines[i] == f"Blurred Image = False\n"
-    i += 1
-    assert lines[i] == f"Mapper = False\n"
-    i += 1
-    assert lines[i] == f"Blurred Mapping Matrix = False\n"
-    i += 1
-    assert lines[i] == f"Inversion Linear Func (Linear Light Profile) Dicts = False\n"
-    i += 1
-    assert lines[i] == f"Curvature Matrix = False\n"
-    i += 1
-    assert lines[i] == f"Curvature Matrix Mapper Diag = False\n"
-    i += 1
-    assert lines[i] == f"Regularization Matrix = False\n"
-    i += 1
-    assert lines[i] == f"Log Det Regularization Matrix Term = False\n"
-    i += 1
-
-    preloads = ag.Preloads(
-        blurred_image=1,
-        w_tilde=1,
-        use_w_tilde=True,
-        image_plane_mesh_grid_pg_list=1,
-        mapper_list=1,
-        operated_mapping_matrix=1,
-        curvature_matrix=1,
-        regularization_matrix=1,
-        log_det_regularization_matrix_term=1,
-    )
-
-    af.formatter.output_list_of_strings_to_file(
-        file=file_preloads, list_of_strings=preloads.info
-    )
-
-    results = open(file_preloads)
-    lines = results.readlines()
-
-    i = 0
-
-    assert lines[i] == f"W Tilde = True\n"
-    i += 1
-    assert lines[i] == f"Use W Tilde = True\n"
-    i += 1
-    assert lines[i] == f"\n"
-    i += 1
-    assert lines[i] == f"Blurred Image = True\n"
-    i += 1
-    assert lines[i] == f"Mapper = True\n"
-    i += 1
-    assert lines[i] == f"Blurred Mapping Matrix = True\n"
-    i += 1
-    assert lines[i] == f"Inversion Linear Func (Linear Light Profile) Dicts = False\n"
-    i += 1
-    assert lines[i] == f"Curvature Matrix = True\n"
-    i += 1
-    assert lines[i] == f"Curvature Matrix Mapper Diag = False\n"
-    i += 1
-    assert lines[i] == f"Regularization Matrix = True\n"
-    i += 1
-    assert lines[i] == f"Log Det Regularization Matrix Term = True\n"
-    i += 1
diff --git a/test_autogalaxy/config/general.yaml b/test_autogalaxy/config/general.yaml
index 2acec5ad..26e205c6 100644
--- a/test_autogalaxy/config/general.yaml
+++ b/test_autogalaxy/config/general.yaml
@@ -1,6 +1,5 @@
 analysis:
   n_cores: 1
-  preload_attempts: 250
 fits:
   flip_for_ds9: true
 grid:
@@ -39,6 +38,4 @@ structures:
   native_binned_only: false           # If True, data structures are only stored in their native and binned format. This is used to reduce memory usage in autocti.
 test:
   check_likelihood_function: false  # if True, when a search is resumed the likelihood of a previous sample is recalculated to ensure it is consistent with the previous run.
-  check_preloads: false
-  exception_override: true
-  preloads_check_threshold: 1.0     # If the figure of merit of a fit with and without preloads is greater than this threshold, the check preload test fails and an exception raised for a model-fit. 
+  exception_override: true
\ No newline at end of file
diff --git a/test_autogalaxy/config/visualize.yaml b/test_autogalaxy/config/visualize.yaml
index 50c40f0c..2b170d7c 100644
--- a/test_autogalaxy/config/visualize.yaml
+++ b/test_autogalaxy/config/visualize.yaml
@@ -8,7 +8,7 @@ include:
     einstein_radius: true
     half_light_radius: true
   include_2d:
-    border: true
+    border: false
     tangential_caustics: false
     radial_caustics: false
     tangential_critical_curves: false
diff --git a/test_autogalaxy/ellipse/test_dataset_interp.py b/test_autogalaxy/ellipse/test_dataset_interp.py
index a2404869..b0e013f3 100644
--- a/test_autogalaxy/ellipse/test_dataset_interp.py
+++ b/test_autogalaxy/ellipse/test_dataset_interp.py
@@ -3,6 +3,32 @@
 import autogalaxy as ag
 
 
+def test__mask_interp():
+    data = ag.Array2D.no_mask(
+        values=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0
+    )
+    noise_map = ag.Array2D.ones(shape_native=(3, 3), pixel_scales=1.0)
+
+    dataset = ag.Imaging(data=data, noise_map=noise_map)
+
+    interp = ag.DatasetInterp(dataset=dataset)
+
+    assert interp.mask_interp((0.5, 0.5)) == pytest.approx(0.0, 1.0e-4)
+
+    mask = ag.Mask2D(
+        mask=[[False, False, False], [False, True, False], [False, False, False]],
+        pixel_scales=1.0,
+    )
+
+    data = data.apply_mask(mask=mask)
+
+    dataset = ag.Imaging(data=data, noise_map=noise_map)
+
+    interp = ag.DatasetInterp(dataset=dataset)
+
+    assert interp.mask_interp((0.5, 0.5)) == pytest.approx(0.25, 1.0e-4)
+
+
 def test__data_interp():
     data = ag.Array2D.no_mask(
         values=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], pixel_scales=1.0
diff --git a/test_autogalaxy/ellipse/test_fit_ellipse.py b/test_autogalaxy/ellipse/test_fit_ellipse.py
index 028ba98c..7a591880 100644
--- a/test_autogalaxy/ellipse/test_fit_ellipse.py
+++ b/test_autogalaxy/ellipse/test_fit_ellipse.py
@@ -4,49 +4,214 @@
 import autogalaxy as ag
 
 
-def test__data_interp(imaging_7x7):
+@pytest.fixture(name="imaging_lh")
+def make_imaging_lh(imaging_7x7):
+    data = ag.Array2D.ones(shape_native=(7, 7), pixel_scales=(1.0, 1.0))
+
+    data[16] = 1.0
+    data[17] = 2.0
+    data[18] = 3.0
+    data[23] = 4.0
+    data[24] = 5.0
+    data[25] = 6.0
+    data[30] = 7.0
+    data[31] = 8.0
+    data[32] = 9.0
+
+    return ag.Imaging(
+        data=data,
+        noise_map=imaging_7x7.noise_map,
+    )
+
+
+@pytest.fixture(name="imaging_lh_masked")
+def make_imaging_lh_masked(imaging_lh):
+    mask = ag.Mask2D(
+        mask=[
+            [True, True, True, True, True, True, True],
+            [True, True, True, True, True, True, True],
+            [True, True, False, False, False, True, True],
+            [True, True, False, True, False, True, True],
+            [True, True, False, False, False, True, True],
+            [True, True, True, True, True, True, True],
+            [True, True, True, True, True, True, True],
+        ],
+        pixel_scales=1.0,
+    )
+
+    return imaging_lh.apply_mask(mask=mask)
+
+
+def test__points_from_major_axis(imaging_lh):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit._points_from_major_axis[1, 0] == pytest.approx(-0.21232, 1.0e-4)
+    assert fit._points_from_major_axis[1, 1] == pytest.approx(0.068987, 1.0e-4)
+
+    assert fit._points_from_major_axis[4, 0] == pytest.approx(0.16366515, 1.0e-4)
+    assert fit._points_from_major_axis[4, 1] == pytest.approx(0.05317803, 1.0e-4)
+
+
+def test___points_from_major_axis__multipole(imaging_lh):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+
+    multipole = ag.EllipseMultipole(m=4, multipole_comps=(0.2, 0.3))
+
+    fit = ag.FitEllipse(
+        dataset=imaging_lh, ellipse=ellipse_0, multipole_list=[multipole]
+    )
+
+    assert fit._points_from_major_axis[1, 0] == pytest.approx(-0.542453, 1.0e-4)
+    assert fit._points_from_major_axis[1, 1] == pytest.approx(-0.038278334, 1.0e-4)
+
+
+def test__mask_interp(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.mask_interp == pytest.approx([False, False, False, False, False], 1.0e-4)
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.mask_interp == pytest.approx([False, True, True, True, True], 1.0e-4)
+
+
+def test__total_points_interp(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.total_points_interp == 5
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.total_points_interp == 1
+
+
+def test__data_interp(imaging_lh, imaging_lh_masked):
     ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
 
-    fit = ag.FitEllipse(dataset=imaging_7x7, ellipse=ellipse_0)
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.data_interp == pytest.approx(
+        [6.0, 2.45584745, 2.42762725, 5.95433876, 8.16218654], 1.0e-4
+    )
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
 
-    assert fit.data_interp[0] == pytest.approx(1.0, 1.0e-4)
-    assert fit.data_interp[1] == pytest.approx(1.0, 1.0e-4)
+    assert fit.data_interp[0] == pytest.approx(6.0, 1.0e-4)
+    assert np.isnan(fit.data_interp[1:5]).all()
 
 
-def test__noise_map_interp(imaging_7x7):
+def test__noise_map_interp(imaging_lh, imaging_lh_masked):
     ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
 
-    fit = ag.FitEllipse(dataset=imaging_7x7, ellipse=ellipse_0)
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.noise_map_interp == pytest.approx([2.0, 2.0, 2.0, 2.0, 2.0], 1.0e-4)
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
 
     assert fit.noise_map_interp[0] == pytest.approx(2.0, 1.0e-4)
-    assert fit.noise_map_interp[1] == pytest.approx(2.0, 1.0e-4)
+    assert np.isnan(fit.noise_map_interp[1:5]).all()
 
 
-def test__log_likelihood(imaging_7x7):
-    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+def test__signal_to_noise_map_interp(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
 
-    fit = ag.FitEllipse(dataset=imaging_7x7, ellipse=ellipse_0)
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
 
-    assert fit.log_likelihood == pytest.approx(0.0, 1.0e-4)
+    assert fit.signal_to_noise_map_interp == pytest.approx(
+        [3.0, 1.22792372, 1.21381362, 2.97716938, 4.08109327], 1.0e-4
+    )
 
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
 
-def test__points_from_major_axis(imaging_7x7):
-    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+    assert fit.signal_to_noise_map_interp[0] == pytest.approx(3.0, 1.0e-4)
+    assert np.isnan(fit.signal_to_noise_map_interp[1:5]).all()
 
-    fit = ag.FitEllipse(dataset=imaging_7x7, ellipse=ellipse_0)
 
-    assert fit._points_from_major_axis[1, 0] == pytest.approx(-0.21232, 1.0e-4)
-    assert fit._points_from_major_axis[1, 1] == pytest.approx(0.068987, 1.0e-4)
+def test__residual_map(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
 
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
 
-def test___points_from_major_axis__multipole(imaging_7x7):
-    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+    assert fit.residual_map == pytest.approx(
+        [1.0, -2.54415255, -2.57237275, 0.95433876, 3.16218654], 1.0e-4
+    )
 
-    multipole = ag.EllipseMultipole(m=4, multipole_comps=(0.2, 0.3))
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
 
-    fit = ag.FitEllipse(
-        dataset=imaging_7x7, ellipse=ellipse_0, multipole_list=[multipole]
+    assert fit.residual_map[0] == pytest.approx(0.0, 1.0e-4)
+    assert np.isnan(fit.noise_map_interp[1:5]).all()
+
+
+def test__normalized_residual_map(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.normalized_residual_map == pytest.approx(
+        [0.5, -1.27207628, -1.28618638, 0.47716938, 1.58109327], 1.0e-4
     )
 
-    assert fit._points_from_major_axis[1, 0] == pytest.approx(-0.542453, 1.0e-4)
-    assert fit._points_from_major_axis[1, 1] == pytest.approx(-0.038278334, 1.0e-4)
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.normalized_residual_map[0] == pytest.approx(0.0, 1.0e-4)
+    assert np.isnan(fit.noise_map_interp[1:5]).all()
+
+
+def test__chi_squared_map(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.chi_squared_map == pytest.approx(
+        [0.25, 1.61817806, 1.65427539, 0.22769062, 2.49985593], 1.0e-4
+    )
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.chi_squared_map[0] == pytest.approx(0.0, 1.0e-4)
+    assert np.isnan(fit.noise_map_interp[1:5]).all()
+
+
+def test__chi_squared(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.chi_squared == pytest.approx(
+        sum([0.25, 1.61817806, 1.65427539, 0.22769062, 2.49985593]), 1.0e-4
+    )
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.chi_squared == pytest.approx(0.0, 1.0e-4)
+
+
+def test__noise_normalization(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.0, 0.0), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.noise_normalization == pytest.approx(16.120857137, 1.0e-4)
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.noise_normalization == pytest.approx(3.224171427, 1.0e-4)
+
+
+def test__log_likelihood(imaging_lh, imaging_lh_masked):
+    ellipse_0 = ag.Ellipse(centre=(0.0, 0.0), ell_comps=(0.5, 0.5), major_axis=1.0)
+
+    fit = ag.FitEllipse(dataset=imaging_lh, ellipse=ellipse_0)
+
+    assert fit.log_likelihood == pytest.approx(-0.16764008373, 1.0e-4)
+
+    fit = ag.FitEllipse(dataset=imaging_lh_masked, ellipse=ellipse_0)
+
+    assert fit.log_likelihood == pytest.approx(0.0, 1.0e-4)