Repository: Jammy2211/PyAutoGalaxy (language: python)

Issue:
Feature/pix refactor
The basis refactor allowed for the following API:

bulge = af.Model(ag.lp_basis.Basis, light_profile_list=gaussians, regularization=ag.reg.Constant)

galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge)

model = af.Collection(galaxies=af.Collection(galaxy=galaxy))
This made all linear objects associated with an Inversion use the following API, in a group fashion (particularly how regularization works).

This PR creates an analogous Pixelization object which means that pixelizations now use a similar API and are handled in an analogous fashion internally in the source code.

Thorough documentation will follow this PR.


Discussion / hints:

Fix patch:
diff --git a/autogalaxy/__init__.py b/autogalaxy/__init__.py
index e3184fa83..8a602864a 100644
--- a/autogalaxy/__init__.py
+++ b/autogalaxy/__init__.py
@@ -5,7 +5,7 @@
 from autoarray.dataset.interferometer import SettingsInterferometer
 from autoarray.instruments import acs
 from autoarray.instruments import euclid
-from autoarray.inversion import pixelizations as pix
+from autoarray.inversion.pixelization import mesh
 from autoarray.inversion import regularization as reg
 from autoarray.inversion.inversion.settings import SettingsInversion
 from autoarray.inversion.inversion.factory import inversion_from as Inversion
@@ -15,8 +15,10 @@
 from autoarray.inversion.inversion.factory import (
     inversion_interferometer_unpacked_from as InversionInterferometer,
 )
-from autoarray.inversion.mappers.factory import mapper_from as Mapper
-from autoarray.inversion.pixelizations.settings import SettingsPixelization
+from autoarray.inversion.pixelization.pixelization import Pixelization
+from autoarray.inversion.pixelization.mappers.mapper_grids import MapperGrids
+from autoarray.inversion.pixelization.mappers.factory import mapper_from as Mapper
+from autoarray.inversion.pixelization.settings import SettingsPixelization
 from autoarray.mask.mask_1d import Mask1D
 from autoarray.mask.mask_2d import Mask2D
 from autoarray.operators.convolver import Convolver
@@ -34,8 +36,9 @@
 from autoarray.structures.grids.iterate_2d import Grid2DIterate
 from autoarray.structures.grids.irregular_2d import Grid2DIrregular
 from autoarray.structures.grids.irregular_2d import Grid2DIrregularUniform
-from autoarray.structures.grids.grid_2d_pixelization import Grid2DRectangular
-from autoarray.structures.grids.grid_2d_pixelization import Grid2DVoronoi
+from autoarray.structures.mesh.rectangular_2d import Mesh2DRectangular
+from autoarray.structures.mesh.voronoi_2d import Mesh2DVoronoi
+from autoarray.structures.mesh.delaunay_2d import Mesh2DDelaunay
 from autoarray.structures.vectors.uniform import VectorYX2D
 from autoarray.structures.vectors.irregular import VectorYX2DIrregular
 from autoarray.layout.region import Region1D
@@ -68,6 +71,7 @@
 from .hyper import hyper_data
 from .analysis.setup import SetupHyper
 from .plane.plane import Plane
+from .plane.to_inversion import AbstractToInversion
 from .plane.to_inversion import PlaneToInversion
 from .profiles.geometry_profiles import EllProfile
 from .profiles import (
@@ -77,7 +81,8 @@
     light_and_mass_profiles as lmp,
     scaling_relations as sr,
 )
-from .profiles.light_profiles.light_profiles_linear import LightProfileLinearObjFunc
+from .profiles.light_profiles import basis as lp_basis
+from .profiles.light_profiles.light_profiles_linear import LightProfileLinearObjFuncList
 from .profiles.light_profiles import light_profiles_init as lp_init
 from .profiles.light_profiles import light_profiles_linear as lp_linear
 from .profiles.light_profiles import light_profiles_operated as lp_operated
diff --git a/autogalaxy/abstract_fit.py b/autogalaxy/abstract_fit.py
index 90e27c0f3..1ac0ec584 100644
--- a/autogalaxy/abstract_fit.py
+++ b/autogalaxy/abstract_fit.py
@@ -31,12 +31,12 @@ def total_mappers(self) -> int:
         """
         The total number of `Mapper` objects used by the inversion in this fit.
 
-        A mapper is created for each galaxy with a pixelization / regularization object.
+        A mapper is created for each galaxy with a pixelization object.
         """
         # TODO : When we add regularization to basis need to change this to reflect mappers.
 
         return len(
-            list(filter(None, self.model_obj.cls_list_from(cls=aa.reg.Regularization)))
+            list(filter(None, self.model_obj.cls_list_from(cls=aa.Pixelization)))
         )
 
     @property
@@ -51,11 +51,7 @@ def perform_inversion(self) -> bool:
         -------
             A bool which is True if an inversion is performed.
         """
-        if self.model_obj.has(cls=aa.pix.Pixelization) or self.model_obj.has(
-            cls=LightProfileLinear
-        ):
-            return True
-        return False
+        return self.model_obj.perform_inversion
 
     @property
     def w_tilde(self) -> Optional[aa.WTildeImaging]:
@@ -89,18 +85,25 @@ def linear_light_profile_intensity_dict(
         This function returns a dictionary which maps every linear light profile instance to its solved for
         `intensity` value in the inversion, so that the intensity value of every light profile can be accessed.
         """
+
         if self.inversion is None:
             return None
 
-        linear_obj_func_list = self.inversion.linear_obj_func_list
+        linear_obj_func_list = self.inversion.cls_list_from(
+            cls=aa.AbstractLinearObjFuncList
+        )
 
         linear_light_profile_intensity_dict = {}
 
         for linear_obj_func in linear_obj_func_list:
 
-            linear_light_profile_intensity_dict[linear_obj_func.light_profile] = float(
-                self.inversion.reconstruction_dict[linear_obj_func]
-            )
+            reconstruction = self.inversion.reconstruction_dict[linear_obj_func]
+
+            for i, light_profile in enumerate(linear_obj_func.light_profile_list):
+
+                linear_light_profile_intensity_dict[light_profile] = float(
+                    reconstruction[i]
+                )
 
         return linear_light_profile_intensity_dict
 
diff --git a/autogalaxy/aggregator/imaging.py b/autogalaxy/aggregator/imaging.py
index 043d9b348..e4d683263 100644
--- a/autogalaxy/aggregator/imaging.py
+++ b/autogalaxy/aggregator/imaging.py
@@ -29,10 +29,10 @@ def _imaging_from(fit: af.Fit, settings_imaging: Optional[aa.SettingsImaging] =
         settings_imaging.relative_accuracy = None
 
     if hasattr(settings_imaging, "sub_size_inversion"):
-        settings_imaging.sub_size_pixelized = settings_imaging.sub_size_inversion
+        settings_imaging.sub_size_pixelization = settings_imaging.sub_size_inversion
 
     if hasattr(settings_imaging, "grid_inversion_class"):
-        settings_imaging.grid_pixelized_class = settings_imaging.grid_inversion_class
+        settings_imaging.grid_pixelization_class = settings_imaging.grid_inversion_class
 
     imaging = aa.Imaging(
         image=data,
diff --git a/autogalaxy/analysis/model_util.py b/autogalaxy/analysis/model_util.py
index 49a2f44c8..0a8ada3da 100644
--- a/autogalaxy/analysis/model_util.py
+++ b/autogalaxy/analysis/model_util.py
@@ -1,12 +1,12 @@
 import logging
 import numpy as np
 from scipy.stats import norm
-from typing import ClassVar, Dict, List
+from typing import ClassVar, Dict, List, Optional
 
 import autofit as af
 import autoarray as aa
 
-from autoarray.inversion.pixelizations.abstract import AbstractPixelization
+from autoarray.inversion.pixelization.mesh.abstract import AbstractMesh
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 
 from autogalaxy.galaxy.galaxy import HyperGalaxy
@@ -34,37 +34,41 @@ def isinstance_or_prior(obj, cls):
     return False
 
 
-def pixelization_from(model: af.Collection) -> AbstractPixelization:
+def mesh_list_from(model: af.Collection) -> Optional[List[AbstractMesh]]:
     """
-    For a model containing one or more galaxies, inspect its attributes and return the `pixelization` of a galaxy
-    provided one galaxy has a pixelization, otherwise it returns none. There cannot be more than one `Pixelization` in
-    a model.
+    For a model containing one or more galaxies, inspect its attributes and return the list of `mesh`'s of each
+    `pixelization` of all galaxies. If no galaxy has pixelization an empty list is returned.
     
     This function expects that the input model is a `Collection` where the first model-component has the
     name `galaxies`, and is itself a `Collection` of `Galaxy` instances. This is the
     standard API for creating a model in PyAutoGalaxy.
 
-    The result of `pixelization_from` is used by the preloading to determine whether certain parts of a
+    The result of `mesh_from` is used by the preloading to determine whether certain parts of a
     calculation can be cached before the non-linear search begins for efficiency.
 
     Parameters
     ----------
-    model : af.Collection
+    model
         Contains the `galaxies` in the model that will be fitted via the non-linear search.
 
     Returns
     -------
-    aa.pix.Pixelization or None:
-        The `Pixelization` of a galaxy, provided one galaxy has a `Pixelization`.
+    The `mesh` of a galaxy, provided one galaxy has a `mesh`.
     """
 
-    for galaxy in model.galaxies:
-        if hasattr(galaxy, "pixelization"):
-            if galaxy.pixelization is not None:
-                if isinstance(galaxy.pixelization, af.Model):
-                    return galaxy.pixelization.cls
-                else:
-                    return galaxy.pixelization
+    instance = model.instance_from_prior_medians()
+
+    mesh_list = []
+
+    for galaxy in instance.galaxies:
+
+        pixelization_list = galaxy.cls_list_from(cls=aa.Pixelization)
+
+        for pixelization in pixelization_list:
+            if pixelization is not None:
+                mesh_list.append(pixelization.mesh)
+
+    return mesh_list
 
 
 def has_pixelization_from(model: af.Collection) -> bool:
@@ -86,76 +90,115 @@ def has_pixelization_from(model: af.Collection) -> bool:
 
     Returns
     -------
-    aa.pix.Pixelization or None:
+    aa.mesh.Mesh or None:
         The `Pixelization` of a galaxy, provided one galaxy has a `Pixelization`.
     """
-    pixelization = pixelization_from(model=model)
+    mesh_list = mesh_list_from(model=model)
 
-    return pixelization is not None
+    return len(mesh_list) > 0
 
 
 def set_upper_limit_of_pixelization_pixels_prior(
-    model: af.Collection, pixels_in_mask: int
+    model: af.Collection,
+    pixels_in_mask: int,
+    lower_limit_no_pixels_below_mask: int = 10,
 ):
     """
-    If the pixelization being fitted in the hyper-model fit is a `VoronoiBrightnessImage` pixelization, this function
-    sets the upper limit of its `pixels` prior to the number of data points in the mask.
+    If the mesh(es) of pixelizations being fitted in the hyper-model fit is a `VoronoiBrightnessImage` pixelization,
+    this function sets the upper limit of its `pixels` prior to the number of data points in the mask.
 
     This ensures the KMeans algorithm does not raise an exception due to having fewer data points than source pixels.
 
     Parameters
     ----------
-    model : Collection
+    model
         The hyper model used by the hyper-fit, which models hyper-components like a `Pixelization` or `HyperGalaxy`'s.
-    result
-        The result of a previous `Analysis` search whose maximum log likelihood model forms the basis of the hyper model.
+    pixels_in_mask
+        The number of pixels in the mask, which are used to set the upper and lower limits of the priors on the
+        number of pixels in the pixelization.
+    lower_limit_no_pixels_below_mask
+        If the prior lower limit on the pixelization's number of pixels is above the number of pixels in the mask,
+        the number of pixels in the mask below which the lower limit is set.
     """
 
-    if hasattr(model, "galaxies"):
+    # TODO : I'm sorry Rich
 
-        for galaxy in model.galaxies:
+    if not hasattr(model, "galaxies"):
+        return
 
-            try:
-                pixelization = getattr(galaxy, "pixelization")
-            except AttributeError:
-                pixelization = None
+    instance = model.instance_from_prior_medians()
 
-            if pixelization is not None:
+    galaxy_key_list_dict = {}
 
-                if hasattr(pixelization, "pixels"):
+    for galaxy in instance.galaxies:
 
-                    if hasattr(pixelization.pixels, "upper_limit"):
+        key_list = []
 
-                        if pixels_in_mask < pixelization.pixels.upper_limit:
+        for key, value in galaxy.__dict__.items():
 
-                            if (
-                                pixelization.cls is aa.pix.DelaunayBrightnessImage
-                                or aa.pix.VoronoiBrightnessImage
-                                or aa.pix.VoronoiNNBrightnessImage
-                            ):
+            if isinstance(value, aa.Pixelization):
+                key_list.append(key)
 
-                                lower_limit = pixelization.pixels.lower_limit
+            elif isinstance(value, af.Model):
+                if type(value.cls) == type(aa.Pixelization):
+                    key_list.append(key)
 
-                                log_str = "MODIFY BEFORE FIT -  A pixelization's pixel UniformPrior upper limit" \
-                                          "was greater than the number of pixels in the mask. It has been " \
-                                          "reduced to the number of pixels in the mask.\,"
+        galaxy_key_list_dict[galaxy] = key_list
 
-                                if lower_limit > pixels_in_mask:
+    for galaxy in model.galaxies:
+
+        if galaxy in galaxy_key_list_dict:
+
+            key_list = galaxy_key_list_dict[galaxy]
+
+            for key in key_list:
+
+                pixelization = getattr(galaxy, key)
+
+                if pixelization is not None:
 
-                                    lower_limit = pixels_in_mask - 1
+                    mesh = pixelization.mesh
 
-                                    logger.info(
-                                        log_str +
-                                        "MODIFY BEFORE FIT - The pixelization's pixel UniformPrior lower_limit was "
-                                        "also above the number of pixels in the mask, and has been reduced"
-                                        "to the number of pixels in the mask minus 1."
+                    if hasattr(mesh, "pixels"):
+
+                        if hasattr(mesh.pixels, "upper_limit"):
+
+                            if pixels_in_mask < mesh.pixels.upper_limit:
+
+                                if (
+                                    mesh.cls is aa.mesh.DelaunayBrightnessImage
+                                    or aa.mesh.VoronoiBrightnessImage
+                                    or aa.mesh.VoronoiNNBrightnessImage
+                                ):
+
+                                    lower_limit = mesh.pixels.lower_limit
+
+                                    log_str = (
+                                        "MODIFY BEFORE FIT -  A pixelization mesh's pixel UniformPrior upper limit"
+                                        "was greater than the number of pixels in the mask. It has been "
+                                        "reduced to the number of pixels in the mask.\,"
                                     )
-                                else:
-                                    logger.info(log_str)
 
-                                pixelization.pixels = af.UniformPrior(
-                                    lower_limit=lower_limit, upper_limit=pixels_in_mask
-                                )
+                                    if lower_limit > pixels_in_mask:
+
+                                        lower_limit = (
+                                            pixels_in_mask
+                                            - lower_limit_no_pixels_below_mask
+                                        )
+
+                                        logger.info(
+                                            log_str
+                                            + "MODIFY BEFORE FIT - The pixelization's mesh's pixel UniformPrior lower_limit was "
+                                            "also above the number of pixels in the mask, and has been reduced"
+                                            "to the number of pixels in the mask minus 10."
+                                        )
+                                    else:
+                                        logger.info(log_str)
+
+                                    mesh.pixels = af.UniformPrior(
+                                        lower_limit=lower_limit,
+                                        upper_limit=pixels_in_mask,
+                                    )
 
 
 def clean_model_of_hyper_images(model):
@@ -275,7 +318,7 @@ def hyper_inversion_model_from(
     if setup_hyper is None:
         return None
 
-    model = result.instance.as_model((AbstractPixelization, AbstractRegularization))
+    model = result.instance.as_model((AbstractMesh, AbstractRegularization))
 
     if not has_pixelization_from(model=model):
         return None
@@ -338,7 +381,7 @@ def hyper_fit_no_noise(
         name=f"{search_previous.paths.name}__hyper_inversion",
         unique_tag=search_previous.paths.unique_tag,
         number_of_cores=search_previous.number_of_cores,
-        **setup_hyper.search_pixelized_dict,
+        **setup_hyper.search_pixelization_dict,
     )
 
     hyper_inversion_result = search.fit(model=hyper_model_inversion, analysis=analysis)
@@ -439,7 +482,7 @@ def hyper_fit(
 
     try:
         set_upper_limit_of_pixelization_pixels_prior(
-            model=hyper_inversion_model, result=result
+            model=hyper_inversion_model, pixels_in_mask=result.mask.pixels_in_mask
         )
     except AttributeError:
         pass
@@ -449,7 +492,7 @@ def hyper_fit(
         name=f"{search_previous.paths.name}__hyper_inversion",
         unique_tag=search_previous.paths.unique_tag,
         number_of_cores=search_previous.number_of_cores,
-        **setup_hyper.search_pixelized_dict,
+        **setup_hyper.search_pixelization_dict,
     )
 
     hyper_inversion_result = search.fit(model=hyper_inversion_model, analysis=analysis)
@@ -487,7 +530,7 @@ def hyper_model_from(
         model components now free parameters.
     """
 
-    model = result.instance.as_model((AbstractPixelization, AbstractRegularization))
+    model = result.instance.as_model((AbstractMesh, AbstractRegularization))
 
     model = clean_model_of_hyper_images(model=model)
 
@@ -576,7 +619,7 @@ def stochastic_model_from(
         model_classes.append(LightProfile)
 
     if include_pixelization:
-        model_classes.append(AbstractPixelization)
+        model_classes.append(AbstractMesh)
 
     if include_regularization:
         model_classes.append(AbstractRegularization)
@@ -610,7 +653,7 @@ def stochastic_model_from(
 def stochastic_fit(
     stochastic_model: af.Collection,
     search_cls: ClassVar[af.NonLinearSearch],
-    search_pixelized_dict: Dict,
+    search_pixelization_dict: Dict,
     result: "ResultDataset",
     analysis: "AnalysisDataset",
     search_previous: af.NonLinearSearch,
@@ -657,7 +700,7 @@ def stochastic_fit(
         name=name,
         unique_tag=search_previous.paths.unique_tag,
         number_of_cores=search_previous.number_of_cores,
-        **search_pixelized_dict,
+        **search_pixelization_dict,
     )
 
     stochastic_result = search.fit(
diff --git a/autogalaxy/analysis/preloads.py b/autogalaxy/analysis/preloads.py
index 531a9d049..37f3d5818 100644
--- a/autogalaxy/analysis/preloads.py
+++ b/autogalaxy/analysis/preloads.py
@@ -130,9 +130,8 @@ def setup_all_via_fits(cls, fit_0, fit_1) -> "Preloads":
         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.plane.to_inversion.mapper_galaxy_dict_from(
-                grid=fit_0.dataset.grid_pixelized
-            )
+
+            preloads.mapper_galaxy_dict = fit_0.plane_to_inversion.mapper_galaxy_dict
 
         preloads.set_operated_mapping_matrix_with_preloads(fit_0=fit_0, fit_1=fit_1)
         preloads.set_curvature_matrix(fit_0=fit_0, fit_1=fit_1)
@@ -199,7 +198,7 @@ def reset_all(self):
         self.sparse_image_plane_grid_pg_list = None
         self.relocated_grid = None
         self.mapper = None
-        self.blurred_mapping_matrix = None
+        self.operated_mapping_matrix = None
         self.curvature_matrix_preload = None
         self.curvature_matrix_counts = None
         self.curvature_matrix = None
diff --git a/autogalaxy/analysis/setup.py b/autogalaxy/analysis/setup.py
index 69fa2517a..a7a180d5f 100644
--- a/autogalaxy/analysis/setup.py
+++ b/autogalaxy/analysis/setup.py
@@ -15,7 +15,7 @@ def __init__(
         search_inversion_cls: Optional[af.NonLinearSearch] = None,
         search_noise_cls: Optional[af.NonLinearSearch] = None,
         search_bc_cls: Optional[af.NonLinearSearch] = None,
-        search_pixelized_dict: Optional[dict] = None,
+        search_pixelization_dict: Optional[dict] = None,
         search_noise_dict: Optional[dict] = None,
         search_bc_dict: Optional[dict] = None,
     ):
@@ -39,7 +39,7 @@ def __init__(
             noise-map's background component in the model.
         search_inversion_cls
             The non-linear search used by every hyper model-fit search.
-        search_pixelized_dict
+        search_pixelization_dict
             The dictionary of search options for the hyper inversion model-fit searches.
         search_noise_dict
             The dictionary of search options for the hyper noise model-fit searches.
@@ -50,7 +50,7 @@ def __init__(
         self.hyper_galaxy_names = None
 
         self.search_inversion_cls = search_inversion_cls or af.DynestyStatic
-        self.search_pixelized_dict = search_pixelized_dict or {
+        self.search_pixelization_dict = search_pixelization_dict or {
             "nlive": 50,
             "sample": "rstagger",
             "dlogz": 10,
diff --git a/autogalaxy/analysis/visualizer.py b/autogalaxy/analysis/visualizer.py
index 84f387798..f81484d03 100644
--- a/autogalaxy/analysis/visualizer.py
+++ b/autogalaxy/analysis/visualizer.py
@@ -372,8 +372,8 @@ def should_plot(name):
             reconstructed_image=should_plot("reconstructed_image")
         )
 
-        inversion_plotter.figures_2d_of_mapper(
-            mapper_index=0,
+        inversion_plotter.figures_2d_of_pixelization(
+            pixelization_index=0,
             reconstructed_image=should_plot("reconstructed_image"),
             reconstruction=should_plot("reconstruction"),
             errors=should_plot("errors"),
@@ -384,7 +384,10 @@ def should_plot(name):
         )
 
         if should_plot("subplot_inversion"):
-            for mapper_index in range(len(inversion.mapper_list)):
+
+            mapper_list = inversion.cls_list_from(cls=aa.AbstractMapper)
+
+            for mapper_index in range(len(mapper_list)):
                 inversion_plotter.subplot_of_mapper(mapper_index=mapper_index)
 
         if not during_analysis:
@@ -393,8 +396,8 @@ def should_plot(name):
 
                 inversion_plotter.figures_2d(reconstructed_image=True)
 
-                inversion_plotter.figures_2d_of_mapper(
-                    mapper_index=0,
+                inversion_plotter.figures_2d_of_pixelization(
+                    pixelization_index=0,
                     reconstructed_image=True,
                     reconstruction=True,
                     errors=True,
diff --git a/autogalaxy/config/general.ini b/autogalaxy/config/general.ini
index 8787c89e5..2a8025ade 100644
--- a/autogalaxy/config/general.ini
+++ b/autogalaxy/config/general.ini
@@ -4,6 +4,7 @@ log_file=output.log
 log_level=INFO
 samples_to_csv=False
 model_results_decimal_places=3
+info_whitespace_length=80
 remove_files=False
 force_pickle_overwrite=False
 identifier_version=4
diff --git a/autogalaxy/config/priors/basis.json b/autogalaxy/config/priors/basis.json
new file mode 100644
index 000000000..52b2562cd
--- /dev/null
+++ b/autogalaxy/config/priors/basis.json
@@ -0,0 +1,1 @@
+{"Basis": {}}
\ No newline at end of file
diff --git a/autogalaxy/config/priors/pixelizations.json b/autogalaxy/config/priors/mesh.json
similarity index 100%
rename from autogalaxy/config/priors/pixelizations.json
rename to autogalaxy/config/priors/mesh.json
diff --git a/autogalaxy/exc.py b/autogalaxy/exc.py
index 7205f5b51..45965eeea 100644
--- a/autogalaxy/exc.py
+++ b/autogalaxy/exc.py
@@ -50,7 +50,7 @@ class PixelizationException(af.exc.FitException):
     """
     Raises exceptions associated with the `inversion/pixelization` modules and `Pixelization` classes.
 
-    For example if a `Rectangular` pixelization has dimensions below 3x3.
+    For example if a `Rectangular` mesh has dimensions below 3x3.
 
     This exception overwrites `autoarray.exc.PixelizationException` in order to add a `FitException`. This means that
     if this exception is raised during a model-fit in the analysis class's `log_likelihood_function` that model
diff --git a/autogalaxy/fixtures.py b/autogalaxy/fixtures.py
index 1720a3a15..d320465b8 100644
--- a/autogalaxy/fixtures.py
+++ b/autogalaxy/fixtures.py
@@ -128,7 +128,7 @@ def make_plane_x2_gal_7x7():
 def make_plane_x2_galaxy_inversion_7x7():
     source_gal_inversion = ag.Galaxy(
         redshift=1.0,
-        pixelization=ag.pix.Rectangular(),
+        pixelization=ag.mesh.Rectangular(),
         regularization=ag.reg.Constant(),
     )
 
@@ -285,6 +285,6 @@ def make_include_2d_all():
         critical_curves=False,
         caustics=False,
         multiple_images=False,
-        mapper_source_pixelization_grid=True,
-        mapper_data_pixelization_grid=True,
+        mapper_source_mesh_grid=True,
+        mapper_data_mesh_grid=True,
     )
diff --git a/autogalaxy/galaxy/galaxy.py b/autogalaxy/galaxy/galaxy.py
index 365048001..b64fcfc78 100644
--- a/autogalaxy/galaxy/galaxy.py
+++ b/autogalaxy/galaxy/galaxy.py
@@ -1,13 +1,10 @@
 from itertools import count
-from typing import Dict, List, Optional, Type, Union
-
 import numpy as np
+from typing import Dict, List, Optional, Type, Union
 
 import autoarray as aa
 import autofit as af
 
-from autoarray.inversion.pixelizations.abstract import AbstractPixelization
-from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoconf.dictable import Dictable
 from autogalaxy import exc
 from autogalaxy.operate.deflections import OperateDeflections
@@ -15,9 +12,6 @@
 from autogalaxy.profiles.geometry_profiles import GeometryProfile
 from autogalaxy.profiles.light_profiles.light_profiles import LightProfile
 from autogalaxy.profiles.light_profiles.light_profiles_linear import LightProfileLinear
-from autogalaxy.profiles.light_profiles.light_profiles_operated import (
-    LightProfileOperated,
-)
 from autogalaxy.profiles.mass_profiles import MassProfile
 
 
@@ -27,12 +21,7 @@ class Galaxy(af.ModelObject, OperateImageList, OperateDeflections, Dictable):
     """
 
     def __init__(
-        self,
-        redshift: float,
-        pixelization: Optional[AbstractPixelization] = None,
-        regularization: Optional[AbstractRegularization] = None,
-        hyper_galaxy: Optional["HyperGalaxy"] = None,
-        **kwargs,
+        self, redshift: float, hyper_galaxy: Optional["HyperGalaxy"] = None, **kwargs
     ):
         """
         Class representing a galaxy, which is composed of attributes used for fitting hyper_galaxies (e.g. light profiles, \
@@ -44,12 +33,8 @@ def __init__(
         ----------
         redshift
             The redshift of the galaxy.
-        pixelization : inversion.Pixelization
+        pixelization
             The pixelization of the galaxy used to reconstruct an observed image using an inversion.
-        regularization : inversion.Regularization
-            The regularization of the pixel-grid used to reconstruct an observed using an inversion.
-        hyper_galaxy
-            The hyper_galaxies-parameters of the hyper_galaxies-galaxy, which is used for performing a hyper_galaxies-analysis on the noise-map.
             
         Attributes
         ----------
@@ -85,18 +70,6 @@ def __init__(
 
             setattr(self, name, val)
 
-        self.pixelization = pixelization
-        self.regularization = regularization
-
-        if pixelization is not None and regularization is None:
-            raise exc.GalaxyException(
-                "If the galaxy has a pixelization, it must also have a regularization."
-            )
-        if pixelization is None and regularization is not None:
-            raise exc.GalaxyException(
-                "If the galaxy has a regularization, it must also have a pixelization."
-            )
-
         self.hyper_galaxy = hyper_galaxy
 
     def __hash__(self):
@@ -104,36 +77,48 @@ def __hash__(self):
 
     def __repr__(self):
         string = "Redshift: {}".format(self.redshift)
-        if self.pixelization:
-            string += "\nPixelization:\n{}".format(str(self.pixelization))
-        if self.regularization:
-            string += "\nRegularization:\n{}".format(str(self.regularization))
-        if self.hyper_galaxy:
-            string += "\nHyper Galaxy:\n{}".format(str(self.hyper_galaxy))
-        if self.cls_list_from(cls=LightProfile):
+
+        if self.has(cls=LightProfile):
             string += "\nLight Profiles:\n{}".format(
                 "\n".join(map(str, self.cls_list_from(cls=LightProfile)))
             )
+
         if self.has(cls=MassProfile):
             string += "\nMass Profiles:\n{}".format(
                 "\n".join(map(str, self.cls_list_from(cls=MassProfile)))
             )
+
+        if self.has(cls=aa.Pixelization):
+            string += "\nPixelization:\n{}".format(str(self.pixelization))
+
+        if self.hyper_galaxy:
+            string += "\nHyper Galaxy:\n{}".format(str(self.hyper_galaxy))
+
         return string
 
     def __eq__(self, other):
         return all(
             (
                 isinstance(other, Galaxy),
-                self.pixelization == other.pixelization,
                 self.redshift == other.redshift,
-                self.hyper_galaxy == other.hyper_galaxy,
                 self.cls_list_from(cls=LightProfile)
                 == other.cls_list_from(cls=LightProfile),
                 self.cls_list_from(cls=MassProfile)
                 == other.cls_list_from(cls=MassProfile),
+                self.cls_list_from(cls=aa.Pixelization)
+                == other.cls_list_from(cls=aa.Pixelization),
+                self.hyper_galaxy == other.hyper_galaxy,
             )
         )
 
+    @property
+    def profile_dict(self) -> Dict:
+        return {
+            key: value
+            for key, value in self.__dict__.items()
+            if isinstance(value, GeometryProfile)
+        }
+
     def dict(self) -> Dict:
         return {
             **{name: profile.dict() for name, profile in self.profile_dict.items()},
@@ -153,17 +138,20 @@ def cls_list_from(self, cls: Type, cls_filtered: Optional[Type] = None) -> List:
         - If `cls=ag.lp.LightProfile` and `cls_filtered=ag.lp.LightProfileLinear`, a list of all light profiles
         excluding those which are linear light profiles will be returned.
 
+        Parameters
+        ----------
+        cls
+            The type of class that a list of instances of this class in the galaxy are returned for.
+        cls_filtered
+            A class type which is filtered and removed from the class list.
+
         Returns
         -------
             The list of objects in the galaxy that inherit from input `cls`.
         """
-        if cls_filtered is not None:
-            return [
-                value
-                for value in self.__dict__.values()
-                if isinstance(value, cls) and not isinstance(value, cls_filtered)
-            ]
-        return [value for value in self.__dict__.values() if isinstance(value, cls)]
+        return aa.util.misc.cls_list_from(
+            values=self.__dict__.values(), cls=cls, cls_filtered=cls_filtered
+        )
 
     def radial_projected_shape_slim_from(self, grid: aa.type.Grid1D2DLike) -> int:
         """
@@ -587,14 +575,6 @@ def mass_angular_within_circle_from(self, radius: float):
                 "You cannot perform a mass-based calculation on a galaxy which does not have a mass-profile"
             )
 
-    @property
-    def profile_dict(self) -> Dict:
-        return {
-            key: value
-            for key, value in self.__dict__.items()
-            if isinstance(value, GeometryProfile)
-        }
-
 
 class HyperGalaxy:
     _ids = count()
diff --git a/autogalaxy/gui/scribbler.py b/autogalaxy/gui/scribbler.py
index 13bd7c0ab..686dbc13a 100644
--- a/autogalaxy/gui/scribbler.py
+++ b/autogalaxy/gui/scribbler.py
@@ -2,6 +2,7 @@
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
+import matplotlib.colors as colors
 
 from skimage.transform import rescale
 
diff --git a/autogalaxy/imaging/fit_imaging.py b/autogalaxy/imaging/fit_imaging.py
index 56db6192e..432432b12 100644
--- a/autogalaxy/imaging/fit_imaging.py
+++ b/autogalaxy/imaging/fit_imaging.py
@@ -147,8 +147,21 @@ def profile_subtracted_image(self) -> aa.Array2D:
         """
         return self.image - self.blurred_image
 
+    @property
+    def plane_to_inversion(self) -> PlaneToInversion:
+        return PlaneToInversion(
+            plane=self.plane,
+            dataset=self.dataset,
+            data=self.profile_subtracted_image,
+            noise_map=self.noise_map,
+            w_tilde=self.w_tilde,
+            settings_pixelization=self.settings_pixelization,
+            settings_inversion=self.settings_inversion,
+            preloads=self.preloads,
+        )
+
     @cached_property
-    def inversion(self) -> Optional[aa.Inversion]:
+    def inversion(self) -> Optional[aa.AbstractInversion]:
         """
         If the plane has linear objects which are used to fit the data (e.g. a linear light profile / pixelization)
         this function returns a linear inversion, where the flux values of these objects (e.g. the `intensity`
@@ -159,18 +172,7 @@ def inversion(self) -> Optional[aa.Inversion]:
         """
 
         if self.perform_inversion:
-
-            plane_to_inversion = PlaneToInversion(plane=self.plane)
-
-            return plane_to_inversion.inversion_imaging_from(
-                dataset=self.dataset,
-                image=self.profile_subtracted_image,
-                noise_map=self.noise_map,
-                w_tilde=self.w_tilde,
-                settings_pixelization=self.settings_pixelization,
-                settings_inversion=self.settings_inversion,
-                preloads=self.preloads,
-            )
+            return self.plane_to_inversion.inversion
 
     @property
     def model_data(self) -> aa.Array2D:
@@ -183,9 +185,7 @@ def model_data(self) -> aa.Array2D:
         If a inversion is included it is the sum of this image and the inversion's reconstruction of the image.
         """
 
-        if self.plane.has(cls=aa.pix.Pixelization) or self.plane.has(
-            cls=LightProfileLinear
-        ):
+        if self.perform_inversion:
 
             return self.blurred_image + self.inversion.mapped_reconstructed_data
 
diff --git a/autogalaxy/imaging/imaging.py b/autogalaxy/imaging/imaging.py
index 1a6ab05e2..8beb48659 100644
--- a/autogalaxy/imaging/imaging.py
+++ b/autogalaxy/imaging/imaging.py
@@ -80,7 +80,8 @@ def via_plane_from(self, plane, grid, name=None):
         )
 
     def via_galaxies_from(self, galaxies, grid, name=None):
-        """Simulate imaging data for this data, as follows:
+        """
+        Simulate imaging data for this data, as follows:
 
         1)  Setup the image-plane grid of the Imaging arrays, which defines the coordinates used for the ray-tracing.
 
@@ -93,7 +94,8 @@ def via_galaxies_from(self, galaxies, grid, name=None):
         4) Plot the image using Matplotlib, if the plot_imaging bool is True.
 
         5) Output the dataset to .fits format if a dataset_path and data_name are specified. Otherwise, return the simulated \
-           imaging data instance."""
+           imaging data instance.
+           """
 
         plane = Plane(
             redshift=float(np.mean([galaxy.redshift for galaxy in galaxies])),
diff --git a/autogalaxy/imaging/plot/fit_imaging_plotters.py b/autogalaxy/imaging/plot/fit_imaging_plotters.py
index e83b09eec..db544ff57 100644
--- a/autogalaxy/imaging/plot/fit_imaging_plotters.py
+++ b/autogalaxy/imaging/plot/fit_imaging_plotters.py
@@ -202,9 +202,9 @@ def subplot_of_galaxies(self, galaxy_index: Optional[int] = None):
             )
             self.figures_2d_of_galaxies(galaxy_index=galaxy_index, model_image=True)
 
-            if self.plane.has(cls=aa.pix.Pixelization):
-                self.inversion_plotter.figures_2d_of_mapper(
-                    mapper_index=0, reconstruction=True
+            if self.plane.has(cls=aa.Pixelization):
+                self.inversion_plotter.figures_2d_of_pixelization(
+                    pixelization_index=0, reconstruction=True
                 )
 
             self.mat_plot_2d.output.subplot_to_figure(
diff --git a/autogalaxy/interferometer/fit_interferometer.py b/autogalaxy/interferometer/fit_interferometer.py
index 5fb10eb36..c32144377 100644
--- a/autogalaxy/interferometer/fit_interferometer.py
+++ b/autogalaxy/interferometer/fit_interferometer.py
@@ -130,8 +130,21 @@ def profile_subtracted_visibilities(self) -> aa.Visibilities:
         """
         return self.visibilities - self.profile_visibilities
 
+    @property
+    def plane_to_inversion(self) -> PlaneToInversion:
+        return PlaneToInversion(
+            plane=self.plane,
+            dataset=self.dataset,
+            data=self.profile_subtracted_visibilities,
+            noise_map=self.noise_map,
+            w_tilde=self.w_tilde,
+            settings_pixelization=self.settings_pixelization,
+            settings_inversion=self.settings_inversion,
+            preloads=self.preloads,
+        )
+
     @cached_property
-    def inversion(self) -> Optional[aa.Inversion]:
+    def inversion(self) -> Optional[aa.AbstractInversion]:
         """
         If the plane has linear objects which are used to fit the data (e.g. a linear light profile / pixelization)
         this function returns a linear inversion, where the flux values of these objects (e.g. the `intensity`
@@ -142,17 +155,7 @@ def inversion(self) -> Optional[aa.Inversion]:
         """
         if self.perform_inversion:
 
-            plane_to_inversion = PlaneToInversion(plane=self.plane)
-
-            return plane_to_inversion.inversion_interferometer_from(
-                dataset=self.dataset,
-                visibilities=self.profile_subtracted_visibilities,
-                noise_map=self.noise_map,
-                w_tilde=self.w_tilde,
-                settings_pixelization=self.settings_pixelization,
-                settings_inversion=self.settings_inversion,
-                preloads=self.preloads,
-            )
+            return self.plane_to_inversion.inversion
 
     @property
     def model_data(self) -> aa.Visibilities:
@@ -165,9 +168,7 @@ def model_data(self) -> aa.Visibilities:
         If a inversion is included it is the sum of these visibilities and the inversion's reconstructed visibilities.
         """
 
-        if self.plane.has(cls=aa.pix.Pixelization) or self.plane.has(
-            cls=LightProfileLinear
-        ):
+        if self.perform_inversion:
 
             return self.profile_visibilities + self.inversion.mapped_reconstructed_data
 
diff --git a/autogalaxy/interferometer/plot/fit_interferometer_plotters.py b/autogalaxy/interferometer/plot/fit_interferometer_plotters.py
index bb06454cf..d343ab546 100644
--- a/autogalaxy/interferometer/plot/fit_interferometer_plotters.py
+++ b/autogalaxy/interferometer/plot/fit_interferometer_plotters.py
@@ -149,11 +149,11 @@ def subplot_fit_real_space(self):
 
             mapper_index = 0
 
-            self.inversion_plotter.figures_2d_of_mapper(
-                mapper_index=mapper_index, reconstructed_image=True
+            self.inversion_plotter.figures_2d_of_pixelization(
+                pixelization_index=mapper_index, reconstructed_image=True
             )
-            self.inversion_plotter.figures_2d_of_mapper(
-                mapper_index=mapper_index, reconstruction=True
+            self.inversion_plotter.figures_2d_of_pixelization(
+                pixelization_index=mapper_index, reconstruction=True
             )
 
             self.mat_plot_2d.output.subplot_to_figure(
diff --git a/autogalaxy/plane/plane.py b/autogalaxy/plane/plane.py
index b742b6f6e..de5bd6aae 100644
--- a/autogalaxy/plane/plane.py
+++ b/autogalaxy/plane/plane.py
@@ -1,6 +1,6 @@
 import json
 import numpy as np
-from typing import Dict, List, Optional, Type
+from typing import Dict, List, Optional, Tuple, Type
 
 import autoarray as aa
 
@@ -9,8 +9,9 @@
 from autogalaxy import exc
 from autogalaxy.galaxy.galaxy import Galaxy
 from autogalaxy.galaxy.galaxy import HyperGalaxy
-from autogalaxy.plane.to_inversion import PlaneToInversion
+from autogalaxy.profiles.light_profiles.basis import Basis
 from autogalaxy.profiles.light_profiles.light_profiles import LightProfile
+from autogalaxy.profiles.light_profiles.light_profiles_linear import LightProfileLinear
 from autogalaxy.profiles.light_profiles.light_profiles_snr import LightProfileSNR
 from autogalaxy.operate.image import OperateImageGalaxies
 from autogalaxy.operate.deflections import OperateDeflections
@@ -70,7 +71,7 @@ def output_to_json(self, file_path: str):
     def galaxy_redshifts(self) -> List[float]:
         return [galaxy.redshift for galaxy in self.galaxies]
 
-    def has(self, cls: Type) -> bool:
+    def has(self, cls: Tuple[Type]) -> bool:
         if self.galaxies is not None:
             return any(list(map(lambda galaxy: galaxy.has(cls=cls), self.galaxies)))
         return False
@@ -97,9 +98,7 @@ def cls_list_from(self, cls: Type) -> List:
         return cls_list
 
     def galaxies_with_cls_list_from(self, cls: Type) -> List[Galaxy]:
-        return list(
-            filter(lambda galaxy: galaxy.has(cls=aa.pix.Pixelization), self.galaxies)
-        )
+        return list(filter(lambda galaxy: galaxy.has(cls=cls), self.galaxies))
 
     @aa.grid_dec.grid_2d_to_structure
     def image_2d_from(
@@ -250,7 +249,7 @@ def traced_grid_from(self, grid: aa.type.Grid2DLike) -> aa.type.Grid2DLike:
     def hyper_galaxies_with_pixelization_image_list(self) -> List[aa.Array2D]:
         return [
             galaxy.hyper_galaxy_image
-            for galaxy in self.galaxies_with_cls_list_from(cls=aa.pix.Pixelization)
+            for galaxy in self.galaxies_with_cls_list_from(cls=aa.Pixelization)
         ]
 
     def hyper_noise_map_from(self, noise_map) -> aa.Array2D:
@@ -322,8 +321,27 @@ def contribution_map_list(self) -> List[aa.Array2D]:
         return contribution_map_list
 
     @property
-    def to_inversion(self):
-        return PlaneToInversion(plane=self)
+    def perform_inversion(self) -> bool:
+        """
+        Returns a bool specifying whether this fit object performs an inversion.
+
+        This is based on whether any of the galaxies in the `model_obj` have a `Pixelization` or `LightProfileLinear`
+        object, in which case an inversion is performed.
+
+        Returns
+        -------
+            A bool which is True if an inversion is performed.
+        """
+        if self.has(cls=(aa.Pixelization, LightProfileLinear)):
+            return True
+        elif self.has(cls=Basis):
+            basis_list = self.cls_list_from(cls=Basis)
+            for basis in basis_list:
+                for light_profile in basis.light_profile_list:
+                    if isinstance(light_profile, LightProfileLinear):
+                        return True
+
+        return False
 
     def extract_attribute(self, cls, attr_name):
         """
diff --git a/autogalaxy/plane/to_inversion.py b/autogalaxy/plane/to_inversion.py
index a575f2abf..17dd2f6c2 100644
--- a/autogalaxy/plane/to_inversion.py
+++ b/autogalaxy/plane/to_inversion.py
@@ -1,120 +1,222 @@
-from typing import Dict, List, Optional, Union
+from typing import Dict, List, Optional, Type, Union
+
+from autoconf import cached_property
 
 import autoarray as aa
 
-from autoarray.inversion.inversion.factory import inversion_imaging_unpacked_from
-from autoarray.inversion.inversion.factory import inversion_interferometer_unpacked_from
+from autoarray.inversion.pixelization.mappers.factory import mapper_from
+from autoarray.inversion.inversion.factory import inversion_unpacked_from
 from autogalaxy.profiles.light_profiles.light_profiles_linear import (
-    LightProfileLinearObjFunc,
+    LightProfileLinearObjFuncList,
 )
+from autogalaxy.profiles.light_profiles.basis import Basis
 from autogalaxy.profiles.light_profiles.light_profiles_linear import LightProfileLinear
 from autogalaxy.galaxy.galaxy import Galaxy
 
 
-class PlaneToInversion:
-    def __init__(self, plane: "Plane"):
+class AbstractToInversion:
+    def cls_light_profile_func_list_galaxy_dict_from(
+        self, cls: Type
+    ) -> Dict[LightProfileLinearObjFuncList, Galaxy]:
+        raise NotImplementedError
 
-        self.plane = plane
+    @cached_property
+    def lp_linear_func_list_galaxy_dict(
+        self,
+    ) -> Dict[LightProfileLinearObjFuncList, Galaxy]:
+
+        raise NotImplementedError
 
-    def lp_linear_func_galaxy_dict_from(
+    @cached_property
+    def mapper_galaxy_dict(self):
+        raise NotImplementedError
+
+    @cached_property
+    def linear_obj_galaxy_dict(
         self,
-        source_grid_slim: aa.type.Grid2DLike,
-        source_blurring_grid_slim: Optional[aa.type.Grid2DLike],
-        convolver: Optional[aa.Convolver] = None,
-    ) -> Dict[LightProfileLinearObjFunc, Galaxy]:
+    ) -> Dict[Union[LightProfileLinearObjFuncList, aa.AbstractMapper], Galaxy]:
+
+        lp_linear_func_galaxy_dict = self.lp_linear_func_list_galaxy_dict
+
+        mapper_galaxy_dict = self.mapper_galaxy_dict
+
+        return {**lp_linear_func_galaxy_dict, **mapper_galaxy_dict}
+
+    @cached_property
+    def linear_obj_list(self) -> List[aa.LinearObj]:
+        return list(self.linear_obj_galaxy_dict.keys())
+
+    @cached_property
+    def regularization_list(self) -> List[aa.AbstractRegularization]:
+        return [linear_obj.regularization for linear_obj in self.linear_obj_list]
+
+
+class PlaneToInversion(AbstractToInversion):
+    def __init__(
+        self,
+        plane: "Plane",
+        dataset: Optional[Union[aa.Imaging, aa.Interferometer]] = None,
+        data: Optional[Union[aa.Array2D, aa.Visibilities]] = None,
+        noise_map: Optional[Union[aa.Array2D, aa.VisibilitiesNoiseMap]] = None,
+        w_tilde: Optional[Union[aa.WTildeImaging, aa.WTildeInterferometer]] = None,
+        grid: Optional[aa.type.Grid2DLike] = None,
+        blurring_grid: Optional[aa.type.Grid2DLike] = None,
+        grid_pixelization: Optional[aa.type.Grid2DLike] = None,
+        settings_pixelization=aa.SettingsPixelization(),
+        settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
+        preloads=aa.Preloads(),
+    ):
 
-        if not self.plane.has(cls=LightProfileLinear):
+        self.plane = plane
+        self.dataset = dataset
+        self.data = data
+        self.noise_map = noise_map
+        self.w_tilde = w_tilde
+
+        if grid is not None:
+            self.grid = grid
+        elif dataset is not None:
+            self.grid = dataset.grid
+        else:
+            self.grid = None
+
+        if blurring_grid is not None:
+            self.blurring_grid = blurring_grid
+        elif dataset is not None:
+            self.blurring_grid = dataset.blurring_grid
+        else:
+            self.blurring_grid = None
+
+        if grid_pixelization is not None:
+            self.grid_pixelization = grid_pixelization
+        elif dataset is not None:
+            self.grid_pixelization = dataset.grid_pixelization
+        else:
+            self.grid_pixelization = None
+
+        self.settings_pixelization = settings_pixelization
+        self.settings_inversion = settings_inversion
+
+        self.preloads = preloads
+
+    def cls_light_profile_func_list_galaxy_dict_from(
+        self, cls: Type
+    ) -> Dict[LightProfileLinearObjFuncList, Galaxy]:
+
+        if not self.plane.has(cls=cls):
             return {}
 
         lp_linear_func_galaxy_dict = {}
 
         for galaxy in self.plane.galaxies:
-            if galaxy.has(cls=LightProfileLinear):
-                for light_profile_linear in galaxy.cls_list_from(
-                    cls=LightProfileLinear
-                ):
+            if galaxy.has(cls=cls):
+                for light_profile in galaxy.cls_list_from(cls=cls):
+
+                    if isinstance(light_profile, LightProfileLinear):
+                        light_profile_list = [light_profile]
+                    else:
+                        light_profile_list = light_profile.light_profile_list
+                        light_profile_list = [
+                            light_profile
+                            for light_profile in light_profile_list
+                            if isinstance(light_profile, LightProfileLinear)
+                        ]
+
+                    if len(light_profile_list) > 0:
+
+                        lp_linear_func = LightProfileLinearObjFuncList(
+                            grid=self.grid,
+                            blurring_grid=self.blurring_grid,
+                            convolver=self.dataset.convolver,
+                            light_profile_list=light_profile_list,
+                            regularization=light_profile.regularization,
+                        )
+
+                        lp_linear_func_galaxy_dict[lp_linear_func] = galaxy
 
-                    lp_linear_func = LightProfileLinearObjFunc(
-                        grid=source_grid_slim,
-                        blurring_grid=source_blurring_grid_slim,
-                        convolver=convolver,
-                        light_profile=light_profile_linear,
-                    )
+        return lp_linear_func_galaxy_dict
 
-                    lp_linear_func_galaxy_dict[lp_linear_func] = galaxy
+    @cached_property
+    def lp_linear_func_list_galaxy_dict(
+        self,
+    ) -> Dict[LightProfileLinearObjFuncList, Galaxy]:
 
-        return lp_linear_func_galaxy_dict
+        lp_linear_light_profile_func_list_galaxy_dict = self.cls_light_profile_func_list_galaxy_dict_from(
+            cls=LightProfileLinear
+        )
 
-    def sparse_image_plane_grid_list_from(
-        self, grid: aa.type.Grid2DLike, settings_pixelization=aa.SettingsPixelization()
-    ) -> Optional[List[aa.type.Grid2DLike]]:
+        lp_basis_func_list_galaxy_dict = self.cls_light_profile_func_list_galaxy_dict_from(
+            cls=Basis
+        )
+
+        return {
+            **lp_linear_light_profile_func_list_galaxy_dict,
+            **lp_basis_func_list_galaxy_dict,
+        }
 
-        if not self.plane.has(cls=aa.pix.Pixelization):
+    @cached_property
+    def sparse_image_plane_grid_list(self,) -> Optional[List[aa.Grid2DSparse]]:
+
+        if not self.plane.has(cls=aa.Pixelization):
             return None
 
         return [
-            pixelization.data_pixelization_grid_from(
-                data_grid_slim=grid,
-                hyper_image=hyper_galaxy_image,
-                settings=settings_pixelization,
+            pixelization.mesh.data_mesh_grid_from(
+                data_grid_slim=self.grid_pixelization,
+                hyper_data=hyper_galaxy_image,
+                settings=self.settings_pixelization,
             )
             for pixelization, hyper_galaxy_image in zip(
-                self.plane.cls_list_from(cls=aa.pix.Pixelization),
+                self.plane.cls_list_from(cls=aa.Pixelization),
                 self.plane.hyper_galaxies_with_pixelization_image_list,
             )
         ]
 
     def mapper_from(
         self,
-        source_grid_slim: aa.type.Grid2DLike,
-        source_pixelization_grid: aa.type.Grid2DLike,
-        pixelization: aa.AbstractPixelization,
+        mesh: aa.AbstractMesh,
+        regularization: aa.AbstractRegularization,
+        source_mesh_grid: aa.Grid2DSparse,
         hyper_galaxy_image: aa.Array2D,
-        data_pixelization_grid: aa.Grid2D = None,
-        settings_pixelization=aa.SettingsPixelization(),
-        preloads=aa.Preloads(),
+        data_mesh_grid: aa.Grid2DSparse = None,
     ) -> aa.AbstractMapper:
 
-        return pixelization.mapper_from(
-            source_grid_slim=source_grid_slim,
-            source_pixelization_grid=source_pixelization_grid,
-            data_pixelization_grid=data_pixelization_grid,
-            hyper_image=hyper_galaxy_image,
-            settings=settings_pixelization,
-            preloads=preloads,
+        mapper_grids = mesh.mapper_grids_from(
+            source_grid_slim=self.grid_pixelization,
+            source_mesh_grid=source_mesh_grid,
+            data_mesh_grid=data_mesh_grid,
+            hyper_data=hyper_galaxy_image,
+            settings=self.settings_pixelization,
+            preloads=self.preloads,
             profiling_dict=self.plane.profiling_dict,
         )
 
-    def mapper_galaxy_dict_from(
-        self,
-        grid: aa.Grid2D,
-        settings_pixelization=aa.SettingsPixelization(),
-        preloads=aa.Preloads(),
-    ) -> Dict[aa.AbstractMapper, Galaxy]:
+        return mapper_from(mapper_grids=mapper_grids, regularization=regularization)
+
+    @cached_property
+    def mapper_galaxy_dict(self) -> Dict[aa.AbstractMapper, Galaxy]:
 
-        if not self.plane.has(cls=aa.pix.Pixelization):
+        if not self.plane.has(cls=aa.Pixelization):
             return {}
 
-        sparse_grid_list = self.sparse_image_plane_grid_list_from(grid=grid)
+        sparse_grid_list = self.sparse_image_plane_grid_list
 
         mapper_galaxy_dict = {}
 
-        pixelization_list = self.plane.cls_list_from(cls=aa.pix.Pixelization)
+        pixelization_list = self.plane.cls_list_from(cls=aa.Pixelization)
         galaxies_with_pixelization_list = self.plane.galaxies_with_cls_list_from(
-            cls=aa.pix.Pixelization
+            cls=aa.Pixelization
         )
         hyper_galaxy_image_list = self.plane.hyper_galaxies_with_pixelization_image_list
 
         for mapper_index in range(len(sparse_grid_list)):
 
             mapper = self.mapper_from(
-                source_grid_slim=grid,
-                source_pixelization_grid=sparse_grid_list[mapper_index],
-                pixelization=pixelization_list[mapper_index],
+                mesh=pixelization_list[mapper_index].mesh,
+                regularization=pixelization_list[mapper_index].regularization,
+                source_mesh_grid=sparse_grid_list[mapper_index],
                 hyper_galaxy_image=hyper_galaxy_image_list[mapper_index],
-                data_pixelization_grid=sparse_grid_list[mapper_index],
-                settings_pixelization=settings_pixelization,
-                preloads=preloads,
+                data_mesh_grid=sparse_grid_list[mapper_index],
             )
 
             galaxy = galaxies_with_pixelization_list[mapper_index]
@@ -123,93 +225,20 @@ def mapper_galaxy_dict_from(
 
         return mapper_galaxy_dict
 
-    def linear_obj_galaxy_dict_from(
-        self,
-        dataset: Union[aa.Imaging, aa.Interferometer],
-        settings_pixelization=aa.SettingsPixelization(),
-        preloads=aa.Preloads(),
-    ) -> Dict[Union[LightProfileLinearObjFunc, aa.AbstractMapper], Galaxy]:
-
-        lp_linear_func_galaxy_dict = self.lp_linear_func_galaxy_dict_from(
-            source_grid_slim=dataset.grid,
-            source_blurring_grid_slim=dataset.blurring_grid,
-            convolver=dataset.convolver,
-        )
-
-        mapper_galaxy_dict = self.mapper_galaxy_dict_from(
-            grid=dataset.grid_pixelized,
-            settings_pixelization=settings_pixelization,
-            preloads=preloads,
-        )
-
-        return {**lp_linear_func_galaxy_dict, **mapper_galaxy_dict}
-
-    def inversion_imaging_from(
-        self,
-        dataset: aa.Imaging,
-        image: aa.Array2D,
-        noise_map: aa.Array2D,
-        w_tilde: aa.WTildeImaging,
-        settings_pixelization: aa.SettingsPixelization = aa.SettingsPixelization(),
-        settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads: aa.Preloads = aa.Preloads(),
-    ):
-
-        linear_obj_galaxy_dict = self.linear_obj_galaxy_dict_from(
-            dataset=dataset,
-            settings_pixelization=settings_pixelization,
-            preloads=preloads,
-        )
-
-        linear_obj_list = list(linear_obj_galaxy_dict.keys())
-
-        inversion = inversion_imaging_unpacked_from(
-            image=image,
-            noise_map=noise_map,
-            convolver=dataset.convolver,
-            w_tilde=w_tilde,
-            linear_obj_list=linear_obj_list,
-            regularization_list=self.plane.cls_list_from(cls=aa.reg.Regularization),
-            settings=settings_inversion,
-            preloads=preloads,
-            profiling_dict=self.plane.profiling_dict,
-        )
-
-        inversion.linear_obj_galaxy_dict = linear_obj_galaxy_dict
-
-        return inversion
-
-    def inversion_interferometer_from(
-        self,
-        dataset: aa.Interferometer,
-        visibilities: aa.Visibilities,
-        noise_map: aa.VisibilitiesNoiseMap,
-        w_tilde,
-        settings_pixelization: aa.SettingsPixelization = aa.SettingsPixelization(),
-        settings_inversion: aa.SettingsInversion = aa.SettingsInversion(),
-        preloads: aa.Preloads = aa.Preloads(),
-    ):
-
-        linear_obj_galaxy_dict = self.linear_obj_galaxy_dict_from(
-            dataset=dataset,
-            settings_pixelization=settings_pixelization,
-            preloads=preloads,
-        )
-
-        linear_obj_list = list(linear_obj_galaxy_dict.keys())
-
-        inversion = inversion_interferometer_unpacked_from(
-            visibilities=visibilities,
-            noise_map=noise_map,
-            transformer=dataset.transformer,
-            w_tilde=w_tilde,
-            linear_obj_list=linear_obj_list,
-            regularization_list=self.plane.cls_list_from(cls=aa.reg.Regularization),
-            settings=settings_inversion,
-            preloads=preloads,
+    @property
+    def inversion(self) -> aa.AbstractInversion:
+
+        inversion = inversion_unpacked_from(
+            dataset=self.dataset,
+            data=self.data,
+            noise_map=self.noise_map,
+            w_tilde=self.w_tilde,
+            linear_obj_list=self.linear_obj_list,
+            settings=self.settings_inversion,
+            preloads=self.preloads,
             profiling_dict=self.plane.profiling_dict,
         )
 
-        inversion.linear_obj_galaxy_dict = linear_obj_galaxy_dict
+        inversion.linear_obj_galaxy_dict = self.linear_obj_galaxy_dict
 
         return inversion
diff --git a/autogalaxy/plot/__init__.py b/autogalaxy/plot/__init__.py
index 26f4b9928..08b9b4a08 100644
--- a/autogalaxy/plot/__init__.py
+++ b/autogalaxy/plot/__init__.py
@@ -36,7 +36,7 @@
     BorderScatter,
     PositionsScatter,
     IndexScatter,
-    PixelizationGridScatter,
+    MeshGridScatter,
     ParallelOverscanPlot,
     SerialPrescanPlot,
     SerialOverscanPlot,
diff --git a/autogalaxy/plot/mat_wrap/include.py b/autogalaxy/plot/mat_wrap/include.py
index e9c07e86d..8ef4bfce9 100644
--- a/autogalaxy/plot/mat_wrap/include.py
+++ b/autogalaxy/plot/mat_wrap/include.py
@@ -33,9 +33,9 @@ def __init__(
         critical_curves=None,
         caustics=None,
         multiple_images=None,
-        mapper_source_pixelization_grid: Optional[bool] = None,
+        mapper_source_mesh_grid: Optional[bool] = None,
         mapper_source_grid_slim: Optional[bool] = None,
-        mapper_data_pixelization_grid=None,
+        mapper_data_mesh_grid=None,
     ):
 
         super().__init__(
@@ -43,9 +43,9 @@ def __init__(
             mask=mask,
             border=border,
             grid=grid,
-            mapper_source_pixelization_grid=mapper_source_pixelization_grid,
+            mapper_source_mesh_grid=mapper_source_mesh_grid,
             mapper_source_grid_slim=mapper_source_grid_slim,
-            mapper_data_pixelization_grid=mapper_data_pixelization_grid,
+            mapper_data_mesh_grid=mapper_data_mesh_grid,
         )
 
         self._positions = positions
diff --git a/autogalaxy/plot/mat_wrap/mat_plot.py b/autogalaxy/plot/mat_wrap/mat_plot.py
index 6857854ab..9e4a20d51 100644
--- a/autogalaxy/plot/mat_wrap/mat_plot.py
+++ b/autogalaxy/plot/mat_wrap/mat_plot.py
@@ -163,7 +163,7 @@ def __init__(
         border_scatter: Optional[aplt.BorderScatter] = None,
         positions_scatter: Optional[aplt.PositionsScatter] = None,
         index_scatter: Optional[aplt.IndexScatter] = None,
-        pixelization_grid_scatter: Optional[aplt.PixelizationGridScatter] = None,
+        mesh_grid_scatter: Optional[aplt.MeshGridScatter] = None,
         light_profile_centres_scatter: Optional[w.LightProfileCentresScatter] = None,
         mass_profile_centres_scatter: Optional[w.MassProfileCentresScatter] = None,
         multiple_images_scatter: Optional[w.MultipleImagesScatter] = None,
@@ -247,7 +247,7 @@ def __init__(
           Scatters specific (y,x) coordinates input as a `Grid2DIrregular` object over the figure.
         index_scatter
           Scatters specific coordinates of an input `Grid2D` based on input values of the `Grid2D`'s 1D or 2D indexes.
-        pixelization_grid_scatter
+        mesh_grid_scatter
           Scatters the `PixelizationGrid` of a `Pixelization` object.
         light_profile_centres_scatter
           Scatters the (y,x) centres of all `LightProfile`'s in the plotted object (e.g. a `Tracer`).
@@ -295,7 +295,7 @@ def __init__(
             grid_scatter=grid_scatter,
             positions_scatter=positions_scatter,
             index_scatter=index_scatter,
-            pixelization_grid_scatter=pixelization_grid_scatter,
+            mesh_grid_scatter=mesh_grid_scatter,
             vector_yx_quiver=vector_yx_quiver,
             patch_overlay=patch_overlay,
             array_overlay=array_overlay,
diff --git a/autogalaxy/plot/mat_wrap/visuals.py b/autogalaxy/plot/mat_wrap/visuals.py
index 3cdbc19ca..e88dbb55e 100644
--- a/autogalaxy/plot/mat_wrap/visuals.py
+++ b/autogalaxy/plot/mat_wrap/visuals.py
@@ -70,7 +70,7 @@ def __init__(
         mask: aa.Mask2D = None,
         positions: Optional[Union[aa.Grid2DIrregular, List[aa.Grid2DIrregular]]] = None,
         grid: Union[aa.Grid2D, aa.Grid2DSparse] = None,
-        pixelization_grid: aa.Grid2D = None,
+        mesh_grid: aa.Grid2D = None,
         vectors: aa.VectorYX2DIrregular = None,
         patches: Union[ptch.Patch] = None,
         array_overlay: aa.Array2D = None,
@@ -92,7 +92,7 @@ def __init__(
             mask=mask,
             positions=positions,
             grid=grid,
-            pixelization_grid=pixelization_grid,
+            mesh_grid=mesh_grid,
             vectors=vectors,
             patches=patches,
             array_overlay=array_overlay,
diff --git a/autogalaxy/profiles/geometry_profiles.py b/autogalaxy/profiles/geometry_profiles.py
index 99a519266..c9eddbfd3 100644
--- a/autogalaxy/profiles/geometry_profiles.py
+++ b/autogalaxy/profiles/geometry_profiles.py
@@ -36,10 +36,7 @@ def has(self, cls: Type) -> bool:
         """
         Does this instance have an attribute which is of type cls?
         """
-        for value in self.__dict__.values():
-            if isinstance(value, cls):
-                return True
-        return False
+        return aa.util.misc.has(values=self.__dict__.values(), cls=cls)
 
     def transform_grid_to_reference_frame(self, grid):
         raise NotImplemented()
diff --git a/autogalaxy/profiles/light_profiles/basis.py b/autogalaxy/profiles/light_profiles/basis.py
new file mode 100644
index 000000000..261ebbcf0
--- /dev/null
+++ b/autogalaxy/profiles/light_profiles/basis.py
@@ -0,0 +1,39 @@
+import numpy as np
+from typing import List, Optional
+
+import autoarray as aa
+
+from autogalaxy.profiles.light_profiles import light_profiles as lp
+from autogalaxy.profiles.light_profiles import light_profiles_linear as lp_linear
+
+
+class Basis(lp.LightProfile):
+    def __init__(
+        self,
+        light_profile_list: List[lp.LightProfile],
+        regularization: Optional[aa.AbstractRegularization] = None,
+    ):
+
+        super().__init__(
+            centre=light_profile_list[0].centre,
+            elliptical_comps=light_profile_list[0].elliptical_comps,
+        )
+
+        self.light_profile_list = light_profile_list
+        self.regularization = regularization
+
+    def image_2d_from(
+        self, grid: aa.type.Grid2DLike, operated_only: Optional[bool] = None
+    ) -> aa.Array2D:
+
+        return sum(self.image_2d_list_from(grid=grid, operated_only=operated_only))
+
+    def image_2d_list_from(
+        self, grid: aa.type.Grid2DLike, operated_only: Optional[bool] = None
+    ) -> List[aa.Array2D]:
+        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],))
+            for light_profile in self.light_profile_list
+        ]
diff --git a/autogalaxy/profiles/light_profiles/light_profiles_linear.py b/autogalaxy/profiles/light_profiles/light_profiles_linear.py
index 86b92881d..2a745650f 100644
--- a/autogalaxy/profiles/light_profiles/light_profiles_linear.py
+++ b/autogalaxy/profiles/light_profiles/light_profiles_linear.py
@@ -1,6 +1,6 @@
 import inspect
 import numpy as np
-from typing import Dict, Optional, Tuple
+from typing import Dict, List, Optional, Tuple
 
 from autoconf import cached_property
 import autoarray as aa
@@ -13,8 +13,14 @@
 from autogalaxy.profiles.light_profiles import light_profiles as lp
 from autogalaxy.profiles import light_and_mass_profiles as lmp
 
+from autogalaxy import exc
+
 
 class LightProfileLinear(lp.LightProfile):
+    @property
+    def regularization(self):
+        return None
+
     @property
     def _intensity(self):
         return 1.0
@@ -87,33 +93,62 @@ def lmp_model_from(self, intensity: float) -> af.Model(lmp.LightMassProfile):
         return af.Model(self.lmp_cls, **parameters_dict)
 
 
-class LightProfileLinearObjFunc(aa.LinearObjFunc):
+class LightProfileLinearObjFuncList(aa.AbstractLinearObjFuncList):
     def __init__(
         self,
         grid: aa.type.Grid1D2DLike,
         blurring_grid: aa.type.Grid1D2DLike,
         convolver: Optional[aa.Convolver],
-        light_profile: LightProfileLinear,
+        light_profile_list: List[LightProfileLinear],
+        regularization=aa.reg.Regularization,
         profiling_dict: Optional[Dict] = None,
     ):
 
-        super().__init__(grid=grid, profiling_dict=profiling_dict)
+        for light_profile in light_profile_list:
+
+            if not isinstance(light_profile, LightProfileLinear):
+
+                raise exc.ProfileException(
+                    """
+                    A light profile that is not a LightProfileLinear object has been input into the
+                    LightProfileLinearObjFuncList object.
+                    
+                    Only children of the LightProfileLinear class can be used in a linear inversion.
+                    """
+                )
+
+        super().__init__(
+            grid=grid, regularization=regularization, profiling_dict=profiling_dict
+        )
 
         self.blurring_grid = blurring_grid
         self.convolver = convolver
-        self.light_profile = light_profile
+        self.light_profile_list = light_profile_list
+
+    @property
+    def pixels(self):
+        return len(self.light_profile_list)
 
     @property
     def mapping_matrix(self) -> np.ndarray:
-        return self.light_profile.image_2d_from(grid=self.grid).binned.slim[:, None]
+
+        mapping_matrix = np.zeros(shape=(self.grid.mask.pixels_in_mask, self.pixels))
+
+        for pixel, light_profile in enumerate(self.light_profile_list):
+
+            image_2d = light_profile.image_2d_from(grid=self.grid).binned.slim
+
+            mapping_matrix[:, pixel] = image_2d
+
+        return mapping_matrix
 
     @cached_property
-    def blurred_mapping_matrix_override(self) -> Optional[np.ndarray]:
+    def operated_mapping_matrix_override(self) -> Optional[np.ndarray]:
         """
         The `LinearEqn` object takes the `mapping_matrix` of each linear object and combines it with the `Convolver`
-        operator to perform a 2D convolution and compute the `blurred_mapping_matrix`.
+        operator to perform a 2D convolution and compute the `operated_mapping_matrix`.
 
-        If this property is overwritten this operation is not performed, with the `blurred_mapping_matrix` output this
+        If this property is overwritten this operation is not performed, with the `operated_mapping_matrix` output this
         property automatically used instead.
 
         This is used for a linear light profile because the in-built mapping matrix convolution does not account for
@@ -127,16 +162,26 @@ def blurred_mapping_matrix_override(self) -> Optional[np.ndarray]:
         performed in the linear equation solvers.
         """
 
-        if isinstance(self.light_profile, LightProfileOperated):
+        if isinstance(self.light_profile_list[0], LightProfileOperated):
             return self.mapping_matrix
 
-        image_2d = self.light_profile.image_2d_from(grid=self.grid)
+        operated_mapping_matrix = np.zeros(
+            shape=(self.grid.mask.pixels_in_mask, self.pixels)
+        )
+
+        for pixel, light_profile in enumerate(self.light_profile_list):
+
+            image_2d = light_profile.image_2d_from(grid=self.grid)
+
+            blurring_image_2d = light_profile.image_2d_from(grid=self.blurring_grid)
+
+            blurred_image_2d = self.convolver.convolve_image(
+                image=image_2d, blurring_image=blurring_image_2d
+            )
 
-        blurring_image_2d = self.light_profile.image_2d_from(grid=self.blurring_grid)
+            operated_mapping_matrix[:, pixel] = blurred_image_2d
 
-        return self.convolver.convolve_image(
-            image=image_2d, blurring_image=blurring_image_2d
-        )[:, None]
+        return operated_mapping_matrix
 
 
 class EllSersic(lp.EllSersic, LightProfileLinear):
diff --git a/autogalaxy/profiles/light_profiles/mock/mock_light_profile.py b/autogalaxy/profiles/light_profiles/mock/mock_light_profile.py
index c6c12a8f2..5a9cb8abf 100644
--- a/autogalaxy/profiles/light_profiles/mock/mock_light_profile.py
+++ b/autogalaxy/profiles/light_profiles/mock/mock_light_profile.py
@@ -1,6 +1,8 @@
 import numpy as np
 from typing import Optional
 
+import autoarray as aa
+
 from autogalaxy.profiles.light_profiles.light_profile_decorators import (
     check_operated_only,
 )
@@ -26,6 +28,7 @@ def __init__(
         self.value = value
         self.value1 = value1
 
+    @aa.grid_dec.grid_2d_to_structure
     @check_operated_only
     def image_2d_from(self, grid, operated_only: Optional[bool] = None):
 
diff --git a/autogalaxy/quantity/dataset_quantity.py b/autogalaxy/quantity/dataset_quantity.py
index bfb9a8956..3f6044e3f 100644
--- a/autogalaxy/quantity/dataset_quantity.py
+++ b/autogalaxy/quantity/dataset_quantity.py
@@ -30,22 +30,22 @@ def __init__(
         grid_class : ag.Grid2D
             The type of grid used to create the image from the `Galaxy` and `Plane`. The options are `Grid2D`,
             and `Grid2DIterate` (see the `Grid2D` documentation for a description of these options).
-        grid_pixelized_class : ag.Grid2D
+        grid_pixelization_class : ag.Grid2D
             The type of grid used to create the grid that maps the `LEq` source pixels to the data's image-pixels.
             The options are `Grid2D` and `Grid2DIterate`.
             (see the `Grid2D` documentation for a description of these options).
         sub_size
-            If the grid and / or grid_pixelized use a `Grid2D`, this sets the sub-size used by the `Grid2D`.
+            If the grid and / or grid_pixelization use a `Grid2D`, this sets the sub-size used by the `Grid2D`.
         fractional_accuracy
-            If the grid and / or grid_pixelized use a `Grid2DIterate`, this sets the fractional accuracy it
+            If the grid and / or grid_pixelization use a `Grid2DIterate`, this sets the fractional accuracy it
             uses when evaluating functions, where the fraction accuracy is the ratio of the values computed using
             two grids at a higher and lower sub-grid size.
         relative_accuracy
-            If the grid and / or grid_pixelized use a `Grid2DIterate`, this sets the relative accuracy it
+            If the grid and / or grid_pixelization use a `Grid2DIterate`, this sets the relative accuracy it
             uses when evaluating functions, where the relative accuracy is the absolute difference of the values
             computed using two grids at a higher and lower sub-grid size.
         sub_steps : [int]
-            If the grid and / or grid_pixelized use a `Grid2DIterate`, this sets the steps the sub-size is increased by
+            If the grid and / or grid_pixelization use a `Grid2DIterate`, this sets the steps the sub-size is increased by
             to meet the fractional accuracy when evaluating functions.
         signal_to_noise_limit
             If input, the dataset's noise-map is rescaled such that no pixel has a signal-to-noise above the
diff --git a/autogalaxy/util/__init__.py b/autogalaxy/util/__init__.py
index a25e2bd1b..9e450c9f8 100644
--- a/autogalaxy/util/__init__.py
+++ b/autogalaxy/util/__init__.py
@@ -7,12 +7,18 @@
 from autoarray.structures.grids import grid_2d_util as grid_2d
 from autoarray.structures.grids import sparse_2d_util as sparse
 from autoarray.fit import fit_util as fit
-from autoarray.inversion.pixelizations import pixelization_util as pixelization
-from autoarray.inversion.mappers import mapper_util as mapper
+from autoarray.inversion.pixelization.mesh import mesh_util as mesh
+from autoarray.inversion.pixelization.mappers import mapper_util as mapper
 from autoarray.inversion.regularization import regularization_util as regularization
-from autoarray.inversion.linear_eqn import leq_util as leq
 from autoarray.inversion.inversion import inversion_util as inversion
+from autoarray.inversion.inversion.imaging import (
+    inversion_imaging_util as inversion_imaging,
+)
+from autoarray.inversion.inversion.interferometer import (
+    inversion_interferometer_util as inversion_interferometer,
+)
 from autoarray.operators import transformer_util as transformer
+from autoarray.util import misc_util as misc
 from autogalaxy.analysis import model_util as model
 
 from autogalaxy.util import error_util as error
diff --git a/docs/api/api.rst b/docs/api/api.rst
index b791c24e5..98da41df5 100644
--- a/docs/api/api.rst
+++ b/docs/api/api.rst
@@ -244,7 +244,7 @@ Plots
     BorderScatter
     PositionsScatter
     IndexScatter
-    PixelizationGridScatter
+    MeshGridScatter
     ParallelOverscanPlot
     SerialPrescanPlot
     SerialOverscanPlot
diff --git a/docs/howtogalaxy/chapter_2_modeling.rst b/docs/howtogalaxy/chapter_2_modeling.rst
index fd0768389..e47f55222 100644
--- a/docs/howtogalaxy/chapter_2_modeling.rst
+++ b/docs/howtogalaxy/chapter_2_modeling.rst
@@ -17,11 +17,14 @@ The chapter contains the following tutorials:
 `Tutorial 4: Dealing with Failure <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_4_dealing_with_failure.ipynb>`_
 - What to do when PyAutoGalaxy finds an inaccurate model.
 
-`Tutorial 5: Masking and Positions <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_5_masking.ipynb>`_
+`Tutorial 5: Linear Profiles <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_5_linear_profiles.ipynb>`_
+- Light profiles which capture complex morphologies in a reduced number of non-linear parameters.
+
+`Tutorial 6: Masking and Positions <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_6_masking.ipynb>`_
 - How to mask your data to improve the model.
 
-`Tutorial 6: Results <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_6_results.ipynb>`_
+`Tutorial 7: Results <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_7_results.ipynb>`_
 - Overview of the results available after successfully fitting a model.
 
-`Tutorial 7: Need for Speed <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_7_need_for_speed.ipynb>`_
+`Tutorial 8: Need for Speed <https://mybinder.org/v2/gh/Jammy2211/autogalaxy_workspace/release?filepath=notebooks/howtogalaxy/chapter_2_modeling/tutorial_8_need_for_speed.ipynb>`_
 - How to fit complex models whilst balancing efficiency and run-time.
diff --git a/docs/index.rst b/docs/index.rst
index 6a6087515..e450d1a40 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -130,7 +130,7 @@ To perform model-fitting, `PyAutoGalaxy` adopts the probabilistic programming
 language `PyAutoFit` (https://github.com/rhayes777/PyAutoFit). `PyAutoFit` allows users to compose a
 model from `LightProfile` and `Galaxy` objects, customize the model parameterization and fit it to data via a
 non-linear search (e.g., `dynesty` [@dynesty], `emcee` [@emcee], `PySwarms` [@pyswarms]). By composing a model with
-`Pixelization` and `Regularization` objects, the galaxy's light is reconstructed using a non-parametric rectangular
+`Pixelization`  objects, the galaxy's light is reconstructed using a non-parametric rectangular
 grid or Voronoi mesh that accounts for irregular galaxy morphologies.
 
 .. code-block:: python
diff --git a/docs/overview/overview_3_modeling.rst b/docs/overview/overview_3_modeling.rst
index 5d443c475..4577d43d5 100644
--- a/docs/overview/overview_3_modeling.rst
+++ b/docs/overview/overview_3_modeling.rst
@@ -362,6 +362,105 @@ model is inferred.
         bar=ag.lp_linear.EllGaussian,
     )
 
+Basis Functions
+---------------
+
+A natural extension of linear light profiles are basis functions, which group many linear light profiles together in
+order to capture complex and irregular structures in a galaxy's emission.
+
+Using a clever model parameterization a basis can be composed which corresponds to just N = 5-10 parameters, making
+model-fitting efficient and robust.
+
+Below, we compose a basis of 10 Gaussians which all share the same `centre` and `elliptical_comps`. Their `sigma`
+values are set via the relation `y = a + (log10(i+1) + b)`, where `i` is the  Gaussian index and `a` and `b` are free
+parameters.
+
+Because `a` and `b` are free parameters (as opposed to `sigma` which can assume many values), we are able to
+compose and fit `Basis` objects which can capture very complex light distributions with just N = 5-10 non-linear
+parameters!
+
+.. code-block:: python
+
+    bulge_a = af.UniformPrior(lower_limit=0.0, upper_limit=0.2)
+    bulge_b = af.UniformPrior(lower_limit=0.0, upper_limit=10.0)
+
+    gaussians_bulge = af.Collection(af.Model(ag.lp_linear.EllGaussian) for _ in range(10))
+
+    for i, gaussian in enumerate(gaussians_bulge):
+
+        gaussian.centre = gaussians_bulge[0].centre
+        gaussian.elliptical_comps = gaussians_bulge[0].elliptical_comps
+        gaussian.sigma = bulge_a + (bulge_b * np.log10(i+1))
+
+    bulge = af.Model(
+        ag.lp_basis.Basis, light_profile_list=gaussians_bulge,
+    )
+
+The bulge's ``info`` attribute describes the basis model composition:
+
+.. code-block:: python
+
+    print(bulge.info)
+
+Below is a snippet of the model, showing that different Gaussians are in the model parameterization:
+
+.. code-block:: bash
+
+    Total Free Parameters = 6
+
+    model                                                                           Basis (N=6)
+        light_profile_list                                                          CollectionPriorModel (N=6)
+            0                                                                       EllGaussian (N=6)
+                sigma                                                               SumPrior (N=2)
+                    other                                                           MultiplePrior (N=1)
+            1                                                                       EllGaussian (N=6)
+                sigma                                                               SumPrior (N=2)
+                    other                                                           MultiplePrior (N=1)
+            2                                                                       EllGaussian (N=6)
+            ...
+            trimmed for conciseness
+            ...
+
+
+    light_profile_list
+        0
+            centre
+                centre_0                                                            GaussianPrior, mean = 0.0, sigma = 0.3
+                centre_1                                                            GaussianPrior, mean = 0.0, sigma = 0.3
+            elliptical_comps
+                elliptical_comps_0                                                  GaussianPrior, mean = 0.0, sigma = 0.3
+                elliptical_comps_1                                                  GaussianPrior, mean = 0.0, sigma = 0.3
+            sigma
+                bulge_a                                                             UniformPrior, lower_limit = 0.0, upper_limit = 0.2
+                other
+                    bulge_b                                                         UniformPrior, lower_limit = 0.0, upper_limit = 10.0
+                    other                                                           0.0
+        1
+            centre
+                centre_0                                                            GaussianPrior, mean = 0.0, sigma = 0.3
+                centre_1                                                            GaussianPrior, mean = 0.0, sigma = 0.3
+            elliptical_comps
+                elliptical_comps_0                                                  GaussianPrior, mean = 0.0, sigma = 0.3
+                elliptical_comps_1                                                  GaussianPrior, mean = 0.0, sigma = 0.3
+            sigma
+                bulge_a                                                             UniformPrior, lower_limit = 0.0, upper_limit = 0.2
+                other
+                    bulge_b                                                         UniformPrior, lower_limit = 0.0, upper_limit = 10.0
+                    other                                                           0.3010299956639812
+        2
+        ...
+        trimmed for conciseness
+        ...
+
+**PyAutoGalaxy** can also apply Bayesian regularization to Basis functions, which smooths the linear light profiles
+(e.g. the Gaussians) in order to prevent over-fitting noise.
+
+.. code-block:: python
+
+    bulge = af.Model(
+        al.lp_basis.Basis, light_profile_list=gaussians_lens, regularization=al.reg.Constant
+    )
+
 Wrap-Up
 -------
 
diff --git a/docs/overview/overview_5_pixelizations.rst b/docs/overview/overview_5_pixelizations.rst
index 647cc1a71..d96553bb1 100644
--- a/docs/overview/overview_5_pixelizations.rst
+++ b/docs/overview/overview_5_pixelizations.rst
@@ -26,8 +26,8 @@ To fit this image with an ``Inversion``, we first mask the ``Imaging`` object:
 
    imaging = imaging.apply_mask(mask=mask_2d)
 
-To reconstruct the galaxy using a pixel-grid, we simply pass it the ``Pixelization`` class we want to reconstruct its
-light using.
+To reconstruct the galaxy on a pixel-grid, called a mesh, we simply pass it the ``Mesh`` class we want to reconstruct its
+light on.
 
 We also pass a ``Regularization`` scheme which applies a smoothness prior on the source reconstruction.
 
@@ -35,13 +35,14 @@ Below, we use a ``Rectangular`` pixelization with resolution 50 x 50 and a ``Con
 
 .. code-block:: python
 
-    galaxy = ag.Galaxy(
-        redshift=1.0,
-        pixelization=ag.pix.Rectangular(shape=(50, 50)),
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(50, 50)),
         regularization=ag.reg.Constant(coefficient=1.0),
     )
 
-Now that our galaxy has a `Pixelization` and `Regularization`, we are able to fit the data using it in the
+    galaxy = ag.Galaxy(redshift=1.0, pixelization=pixelization)
+
+Now that our galaxy has a `Pixelization`, we are able to fit the data using it in the
 same way as before, by simply passing the galaxy to a `Plane` and using this `Plane` to create a `FitImaging`
 object.
 
diff --git a/paper/paper.md b/paper/paper.md
index b6b657b5e..edffbf7d5 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -133,7 +133,7 @@ To perform model-fitting, `PyAutoGalaxy` adopts the probabilistic programming
 language `PyAutoFit` (https://github.com/rhayes777/PyAutoFit). `PyAutoFit` allows users to compose a 
 model from `LightProfile` and `Galaxy` objects, customize the model parameterization and fit it to data via a 
 non-linear search, for example `dynesty` [@dynesty], `emcee` [@emcee] or `PySwarms` [@pyswarms]. By composing a model with 
-`Pixelization` and `Regularization` objects, the galaxy's light is reconstructed using a non parametric rectangular 
+`Pixelization`  objects, the galaxy's light is reconstructed using a non parametric rectangular 
 grid that accounts for irregular galaxy morphologies. Multiple images of the same galaxy can be fitted simultaneously,
 using models whose parameters vary across wavelength.
 


Test patch:
diff --git a/test_autogalaxy/aggregator/test_aggregator_imaging.py b/test_autogalaxy/aggregator/test_aggregator_imaging.py
index 2d72d8015..b4f5bbb5a 100644
--- a/test_autogalaxy/aggregator/test_aggregator_imaging.py
+++ b/test_autogalaxy/aggregator/test_aggregator_imaging.py
@@ -71,7 +71,7 @@ def test__imaging_generator_from_aggregator(imaging_7x7, mask_2d_7x7, samples, m
     masked_imaging_7x7 = masked_imaging_7x7.apply_settings(
         settings=ag.SettingsImaging(
             grid_class=ag.Grid2DIterate,
-            grid_pixelized_class=ag.Grid2DIterate,
+            grid_pixelization_class=ag.Grid2DIterate,
             fractional_accuracy=0.5,
             sub_steps=[2],
         )
@@ -93,7 +93,7 @@ def test__imaging_generator_from_aggregator(imaging_7x7, mask_2d_7x7, samples, m
     for imaging in imaging_gen:
         assert (imaging.image == masked_imaging_7x7.image).all()
         assert isinstance(imaging.grid, ag.Grid2DIterate)
-        assert isinstance(imaging.grid_pixelized, ag.Grid2DIterate)
+        assert isinstance(imaging.grid_pixelization, ag.Grid2DIterate)
         assert imaging.grid.sub_steps == [2]
         assert imaging.grid.fractional_accuracy == 0.5
 
diff --git a/test_autogalaxy/aggregator/test_aggregator_interferometer.py b/test_autogalaxy/aggregator/test_aggregator_interferometer.py
index 8ed8313ef..906948cb4 100644
--- a/test_autogalaxy/aggregator/test_aggregator_interferometer.py
+++ b/test_autogalaxy/aggregator/test_aggregator_interferometer.py
@@ -81,7 +81,7 @@ def test__interferometer_generator_from_aggregator(
         real_space_mask=mask_2d_7x7,
         settings=ag.SettingsInterferometer(
             grid_class=ag.Grid2DIterate,
-            grid_pixelized_class=ag.Grid2DIterate,
+            grid_pixelization_class=ag.Grid2DIterate,
             fractional_accuracy=0.5,
             sub_steps=[2],
             transformer_class=ag.TransformerDFT,
@@ -105,7 +105,7 @@ def test__interferometer_generator_from_aggregator(
         assert (interferometer.visibilities == interferometer_7.visibilities).all()
         assert (interferometer.real_space_mask == mask_2d_7x7).all()
         assert isinstance(interferometer.grid, ag.Grid2DIterate)
-        assert isinstance(interferometer.grid_pixelized, ag.Grid2DIterate)
+        assert isinstance(interferometer.grid_pixelization, ag.Grid2DIterate)
         assert interferometer.grid.sub_steps == [2]
         assert interferometer.grid.fractional_accuracy == 0.5
         assert isinstance(interferometer.transformer, ag.TransformerDFT)
diff --git a/test_autogalaxy/analysis/test_analysis.py b/test_autogalaxy/analysis/test_analysis.py
index 068b969ee..2ba25e636 100644
--- a/test_autogalaxy/analysis/test_analysis.py
+++ b/test_autogalaxy/analysis/test_analysis.py
@@ -76,15 +76,12 @@ def test__modify_before_fit__kmeans_pixelization_upper_limit_ajusted_based_on_ma
     masked_imaging_7x7
 ):
 
-    pixelization = af.Model(ag.pix.DelaunayBrightnessImage)
+    mesh = af.Model(ag.mesh.DelaunayBrightnessImage)
+    mesh.pixels = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)
 
-    pixelization.pixels = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)
+    pixelization = af.Model(ag.Pixelization, mesh=mesh)
 
-    galaxies = af.Collection(
-        source=ag.Galaxy(
-            redshift=0.5, pixelization=pixelization, regularization=ag.reg.Constant
-        )
-    )
+    galaxies = af.Collection(source=ag.Galaxy(redshift=0.5, pixelization=pixelization))
 
     model = af.Collection(galaxies=galaxies)
 
@@ -92,6 +89,6 @@ def test__modify_before_fit__kmeans_pixelization_upper_limit_ajusted_based_on_ma
 
     analysis.modify_before_fit(paths=af.DirectoryPaths(), model=model)
 
-    assert model.galaxies.source.pixelization.pixels.upper_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.upper_limit == pytest.approx(
         9, 1.0e-4
     )
diff --git a/test_autogalaxy/analysis/test_model_util.py b/test_autogalaxy/analysis/test_model_util.py
index 8842bb702..959ada0b4 100644
--- a/test_autogalaxy/analysis/test_model_util.py
+++ b/test_autogalaxy/analysis/test_model_util.py
@@ -4,87 +4,57 @@
 import autogalaxy as ag
 
 
-def test__pixelization_from_model():
+def test__mesh_list_from_model():
 
     galaxies = af.Collection(galaxy=af.Model(ag.Galaxy, redshift=0.5))
 
-    pixelization = ag.util.model.pixelization_from(
-        model=af.Collection(galaxies=galaxies)
-    )
+    mesh_list = ag.util.model.mesh_list_from(model=af.Collection(galaxies=galaxies))
+
+    assert mesh_list == []
+
+    pixelization_0 = ag.Pixelization(mesh=ag.mesh.Rectangular())
 
-    assert pixelization is None
+    pixelization_1 = ag.Pixelization(mesh=ag.mesh.VoronoiMagnification())
 
     galaxies = af.Collection(
         galaxy=ag.Galaxy(
-            redshift=0.5,
-            pixelization=ag.pix.Rectangular(),
-            regularization=ag.reg.Constant(),
+            redshift=0.5, pixelization_0=pixelization_0, pixelization_1=pixelization_1
         )
     )
 
-    pixelization = ag.util.model.pixelization_from(
-        model=af.Collection(galaxies=galaxies)
-    )
+    mesh_list = ag.util.model.mesh_list_from(model=af.Collection(galaxies=galaxies))
+
+    assert isinstance(mesh_list[0], ag.mesh.Rectangular)
+    assert isinstance(mesh_list[1], ag.mesh.VoronoiMagnification)
 
-    assert isinstance(pixelization, ag.pix.Rectangular)
+    pixelization_0 = af.Model(ag.Pixelization, mesh=ag.mesh.Rectangular)
+
+    pixelization_1 = af.Model(ag.Pixelization, mesh=ag.mesh.VoronoiMagnification)
 
     galaxies = af.Collection(
         galaxy=af.Model(
             ag.Galaxy,
             redshift=0.5,
-            pixelization=ag.pix.Rectangular,
-            regularization=ag.reg.Constant,
+            pixelization_0=pixelization_0,
+            pixelization_1=pixelization_1,
         )
     )
 
-    pixelization = ag.util.model.pixelization_from(
-        model=af.Collection(galaxies=galaxies)
-    )
-
-    assert type(pixelization) == type(ag.pix.Rectangular)
-
-
-# def test__has_pixelization():
-#     galaxies = af.Collection(galaxy=af.Model(ag.Galaxy, redshift=0.5))
-#
-#     assert galaxies.has_model(cls=ag.pix.Pixelization) is False
-#
-#     galaxies = af.Collection(
-#         galaxy=af.Model(
-#             ag.Galaxy,
-#             redshift=0.5,
-#             pixelization=ag.pix.Rectangular(),
-#             regularization=ag.reg.Constant(),
-#         )
-#     )
-#
-#     assert galaxies.has_instance(cls=ag.pix.Pixelization) is True
-#     assert galaxies.has_model(cls=ag.pix.Pixelization) is False
-#
-#     galaxies = af.Collection(
-#         galaxy=af.Model(
-#             ag.Galaxy,
-#             redshift=0.5,
-#             pixelization=ag.pix.Rectangular,
-#             regularization=ag.reg.Constant,
-#         )
-#     )
-#
-#     assert galaxies.has_instance(cls=ag.pix.Pixelization) is False
-#     assert galaxies.has_model(cls=ag.pix.Pixelization) is True
+    mesh_list = ag.util.model.mesh_list_from(model=af.Collection(galaxies=galaxies))
+
+    assert isinstance(mesh_list[0], ag.mesh.Rectangular)
+    assert isinstance(mesh_list[1], ag.mesh.VoronoiMagnification)
 
 
 def test__set_upper_limit_of_pixelization_pixels_prior():
 
-    pixelization = af.Model(ag.pix.DelaunayBrightnessImage)
+    mesh = af.Model(ag.mesh.DelaunayBrightnessImage)
 
-    pixelization.pixels = af.UniformPrior(lower_limit=5.0, upper_limit=10.0)
+    mesh.pixels = af.UniformPrior(lower_limit=5.0, upper_limit=10.0)
 
-    galaxies = af.Collection(
-        source=ag.Galaxy(
-            redshift=0.5, pixelization=pixelization, regularization=ag.reg.Constant
-        )
-    )
+    pixelization = ag.Pixelization(mesh=mesh)
+
+    galaxies = af.Collection(source=ag.Galaxy(redshift=0.5, pixelization=pixelization))
 
     model = af.Collection(galaxies=galaxies)
 
@@ -92,21 +62,31 @@ def test__set_upper_limit_of_pixelization_pixels_prior():
         model=model, pixels_in_mask=12
     )
 
-    assert model.galaxies.source.pixelization.pixels.lower_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.lower_limit == pytest.approx(
         5, 1.0e-4
     )
-    assert model.galaxies.source.pixelization.pixels.upper_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.upper_limit == pytest.approx(
         10, 1.0e-4
     )
 
+    mesh = af.Model(ag.mesh.DelaunayBrightnessImage)
+    mesh.pixels = af.UniformPrior(lower_limit=5.0, upper_limit=10.0)
+    pixelization = af.Model(ag.Pixelization, mesh=mesh)
+    galaxies = af.Collection(source=ag.Galaxy(redshift=0.5, pixelization=pixelization))
+    model = af.Collection(galaxies=galaxies)
+
+    ag.util.model.set_upper_limit_of_pixelization_pixels_prior(
+        model=model, pixels_in_mask=12
+    )
+
     ag.util.model.set_upper_limit_of_pixelization_pixels_prior(
         model=model, pixels_in_mask=8
     )
 
-    assert model.galaxies.source.pixelization.pixels.lower_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.lower_limit == pytest.approx(
         5, 1.0e-4
     )
-    assert model.galaxies.source.pixelization.pixels.upper_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.upper_limit == pytest.approx(
         8, 1.0e-4
     )
 
@@ -114,20 +94,24 @@ def test__set_upper_limit_of_pixelization_pixels_prior():
         model=model, pixels_in_mask=3
     )
 
-    assert model.galaxies.source.pixelization.pixels.lower_limit == pytest.approx(
-        2, 1.0e-4
+    assert model.galaxies.source.pixelization.mesh.pixels.lower_limit == pytest.approx(
+        -7, 1.0e-4
     )
-    assert model.galaxies.source.pixelization.pixels.upper_limit == pytest.approx(
+    assert model.galaxies.source.pixelization.mesh.pixels.upper_limit == pytest.approx(
         3, 1.0e-4
     )
 
-    pixelization = af.Model(ag.pix.VoronoiBrightnessImage)
+    mesh_0 = af.Model(ag.mesh.VoronoiBrightnessImage)
+    mesh_0.pixels = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)
+    pixelization_0 = ag.Pixelization(mesh=mesh_0)
 
-    pixelization.pixels = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)
+    mesh_1 = af.Model(ag.mesh.DelaunayBrightnessImage)
+    mesh_1.pixels = af.UniformPrior(lower_limit=0.0, upper_limit=4.0)
+    pixelization_1 = ag.Pixelization(mesh=mesh_1)
 
     galaxies = af.Collection(
         source=ag.Galaxy(
-            redshift=0.5, pixelization=pixelization, regularization=ag.reg.Constant
+            redshift=0.5, pixelization_0=pixelization_0, pixelization_1=pixelization_1
         )
     )
 
@@ -137,20 +121,26 @@ def test__set_upper_limit_of_pixelization_pixels_prior():
         model=model, pixels_in_mask=6
     )
 
-    assert model.galaxies.source.pixelization.pixels.upper_limit == pytest.approx(
-        6, 1.0e-4
+    assert (
+        model.galaxies.source.pixelization_0.mesh.pixels.upper_limit
+        == pytest.approx(6, 1.0e-4)
+    )
+
+    assert (
+        model.galaxies.source.pixelization_1.mesh.pixels.upper_limit
+        == pytest.approx(4, 1.0e-4)
     )
 
 
 def test__hyper_model_noise_from():
+
+    pixelization = af.Model(
+        ag.Pixelization, mesh=ag.mesh.Rectangular, regularization=ag.reg.Constant
+    )
+
     model = af.Collection(
         galaxies=af.Collection(
-            galaxy=af.Model(
-                ag.Galaxy,
-                redshift=0.5,
-                pixelization=ag.pix.Rectangular,
-                regularization=ag.reg.Constant,
-            ),
+            galaxy=af.Model(ag.Galaxy, redshift=0.5, pixelization=pixelization),
             galaxy_1=af.Model(ag.Galaxy, redshift=1.0, bulge=ag.lp.EllSersic),
         )
     )
@@ -178,11 +168,11 @@ def test__hyper_model_noise_from():
         include_hyper_image_sky=True,
     )
 
-    assert model.galaxies.galaxy.pixelization.cls is ag.pix.Rectangular
-    assert model.galaxies.galaxy.regularization.cls is ag.reg.Constant
+    assert model.galaxies.galaxy.pixelization.mesh.cls is ag.mesh.Rectangular
+    assert model.galaxies.galaxy.pixelization.regularization.cls is ag.reg.Constant
 
     assert model.galaxies.galaxy.pixelization.prior_count == 0
-    assert model.galaxies.galaxy.regularization.prior_count == 0
+    assert model.galaxies.galaxy.pixelization.regularization.prior_count == 0
 
     assert model.galaxies.galaxy_1.bulge.intensity == pytest.approx(1.0, 1.0e-4)
 
@@ -265,14 +255,11 @@ def test__hyper_model_noise_from__adds_hyper_galaxies():
 
 def test__hyper_model_inversion_from():
 
+    pixelization = af.Model(ag.Pixelization, mesh=ag.mesh.Rectangular)
+
     model = af.Collection(
         galaxies=af.Collection(
-            galaxy=af.Model(
-                ag.Galaxy,
-                redshift=0.5,
-                pixelization=ag.pix.Rectangular,
-                regularization=ag.reg.Constant,
-            ),
+            galaxy=af.Model(ag.Galaxy, redshift=0.5, pixelization=pixelization),
             galaxy_1=af.Model(ag.Galaxy, redshift=1.0, bulge=ag.lp.EllSersic),
         )
     )
@@ -285,11 +272,9 @@ def test__hyper_model_inversion_from():
         setup_hyper=ag.SetupHyper(), result=result
     )
 
-    assert isinstance(model.galaxies.galaxy.pixelization, af.Model)
-    assert isinstance(model.galaxies.galaxy.regularization, af.Model)
+    assert isinstance(model.galaxies.galaxy.pixelization.mesh, af.Model)
 
-    assert model.galaxies.galaxy.pixelization.cls is ag.pix.Rectangular
-    assert model.galaxies.galaxy.regularization.cls is ag.reg.Constant
+    assert model.galaxies.galaxy.pixelization.mesh.cls is ag.mesh.Rectangular
     assert model.galaxies.galaxy_1.bulge.intensity == pytest.approx(1.0, 1.0e-4)
 
     assert model.hyper_image_sky is None
@@ -308,8 +293,7 @@ def test__hyper_model_inversion_from():
         include_hyper_image_sky=True,
     )
 
-    assert isinstance(model.galaxies.galaxy.pixelization, af.Model)
-    assert isinstance(model.galaxies.galaxy.regularization, af.Model)
+    assert isinstance(model.galaxies.galaxy.pixelization.mesh, af.Model)
     assert isinstance(model.hyper_image_sky, af.Model)
 
     assert model.hyper_background_noise is None
@@ -335,6 +319,9 @@ def test__hyper_model_inversion_from():
 
 
 def test__hyper_model_inversion_from__adds_hyper_galaxies():
+
+    pixelization = af.Model(ag.Pixelization, mesh=ag.mesh.Rectangular)
+
     model = af.Collection(
         galaxies=af.Collection(
             galaxy_0=af.Model(ag.Galaxy, redshift=0.5),
@@ -342,8 +329,7 @@ def test__hyper_model_inversion_from__adds_hyper_galaxies():
                 ag.Galaxy,
                 redshift=1.0,
                 bulge=ag.lp.EllSersic,
-                pixelization=ag.pix.Rectangular,
-                regularization=ag.reg.Constant,
+                pixelization=pixelization,
             ),
         )
     )
@@ -404,6 +390,13 @@ def test__hyper_model_inversion_from__adds_hyper_galaxies():
 
 
 def test__stochastic_model_from():
+
+    pixelization = af.Model(
+        ag.Pixelization,
+        mesh=ag.mesh.VoronoiBrightnessImage(),
+        regularization=ag.reg.AdaptiveBrightness(),
+    )
+
     model = af.Collection(
         galaxies=af.Collection(
             lens=af.Model(
@@ -412,12 +405,7 @@ def test__stochastic_model_from():
                 light=ag.lp.SphSersic(),
                 mass=ag.mp.SphIsothermal(),
             ),
-            source=af.Model(
-                ag.Galaxy,
-                redshift=1.0,
-                pixelization=ag.pix.VoronoiBrightnessImage(),
-                regularization=ag.reg.AdaptiveBrightness(),
-            ),
+            source=af.Model(ag.Galaxy, redshift=1.0, pixelization=pixelization),
         )
     )
 
@@ -433,15 +421,19 @@ def test__stochastic_model_from():
 
     assert isinstance(model.galaxies.lens.mass.centre, af.TuplePrior)
     assert isinstance(model.galaxies.lens.light.intensity, float)
-    assert isinstance(model.galaxies.source.pixelization.pixels, int)
-    assert isinstance(model.galaxies.source.regularization.inner_coefficient, float)
+    assert isinstance(model.galaxies.source.pixelization.mesh.pixels, int)
+    assert isinstance(
+        model.galaxies.source.pixelization.regularization.inner_coefficient, float
+    )
 
     model = ag.util.model.stochastic_model_from(result=result, include_lens_light=True)
 
     assert isinstance(model.galaxies.lens.mass.centre, af.TuplePrior)
     assert isinstance(model.galaxies.lens.light.intensity, af.LogUniformPrior)
-    assert isinstance(model.galaxies.source.pixelization.pixels, int)
-    assert isinstance(model.galaxies.source.regularization.inner_coefficient, float)
+    assert isinstance(model.galaxies.source.pixelization.mesh.pixels, int)
+    assert isinstance(
+        model.galaxies.source.pixelization.regularization.inner_coefficient, float
+    )
 
     model = ag.util.model.stochastic_model_from(
         result=result, include_pixelization=True
@@ -449,9 +441,10 @@ def test__stochastic_model_from():
 
     assert isinstance(model.galaxies.lens.mass.centre, af.TuplePrior)
     assert isinstance(model.galaxies.lens.light.intensity, float)
-    assert isinstance(model.galaxies.source.pixelization.pixels, af.UniformPrior)
+    assert isinstance(model.galaxies.source.pixelization.mesh.pixels, af.UniformPrior)
     assert not isinstance(
-        model.galaxies.source.regularization.inner_coefficient, af.UniformPrior
+        model.galaxies.source.pixelization.regularization.inner_coefficient,
+        af.UniformPrior,
     )
 
     model = ag.util.model.stochastic_model_from(
@@ -460,7 +453,8 @@ def test__stochastic_model_from():
 
     assert isinstance(model.galaxies.lens.mass.centre, af.TuplePrior)
     assert isinstance(model.galaxies.lens.light.intensity, float)
-    assert isinstance(model.galaxies.source.pixelization.pixels, int)
+    assert isinstance(model.galaxies.source.pixelization.mesh.pixels, int)
     assert isinstance(
-        model.galaxies.source.regularization.inner_coefficient, af.UniformPrior
+        model.galaxies.source.pixelization.regularization.inner_coefficient,
+        af.UniformPrior,
     )
diff --git a/test_autogalaxy/analysis/test_result.py b/test_autogalaxy/analysis/test_result.py
index 6b0791ae9..2774df84f 100644
--- a/test_autogalaxy/analysis/test_result.py
+++ b/test_autogalaxy/analysis/test_result.py
@@ -45,11 +45,10 @@ def test__max_log_likelihood_plane_available_as_result(analysis_imaging_7x7):
 
 
 def test__results_include_pixelization__available_as_property(analysis_imaging_7x7):
-    source = ag.Galaxy(
-        redshift=1.0,
-        pixelization=ag.m.MockPixelization(mapper=1),
-        regularization=ag.m.MockRegularization(),
-    )
+
+    pixelization = ag.m.MockPixelization(mapper=1)
+
+    source = ag.Galaxy(redshift=1.0, pixelization=pixelization)
 
     max_log_likelihood_plane = ag.Plane(galaxies=[source])
 
@@ -60,6 +59,6 @@ def test__results_include_pixelization__available_as_property(analysis_imaging_7
     )
 
     assert isinstance(
-        result.cls_list_from(cls=ag.pix.Pixelization)[0], ag.m.MockPixelization
+        result.cls_list_from(cls=ag.Pixelization)[0], ag.m.MockPixelization
     )
-    assert result.cls_list_from(cls=ag.pix.Pixelization)[0].mapper == 1
+    assert result.cls_list_from(cls=ag.Pixelization)[0].mapper == 1
diff --git a/test_autogalaxy/analysis/test_setup.py b/test_autogalaxy/analysis/test_setup.py
index cd063d42e..91e8d6666 100644
--- a/test_autogalaxy/analysis/test_setup.py
+++ b/test_autogalaxy/analysis/test_setup.py
@@ -5,9 +5,9 @@
 class TestSetupHyper:
     def test__hyper_search(self):
 
-        setup = ag.SetupHyper(search_inversion_cls=None, search_pixelized_dict=None)
+        setup = ag.SetupHyper(search_inversion_cls=None, search_pixelization_dict=None)
         assert setup.search_inversion_cls == af.DynestyStatic
-        assert setup.search_pixelized_dict == {
+        assert setup.search_pixelization_dict == {
             "nlive": 50,
             "sample": "rstagger",
             "dlogz": 10,
@@ -15,10 +15,10 @@ def test__hyper_search(self):
 
         setup = ag.SetupHyper(
             search_inversion_cls=af.DynestyDynamic,
-            search_pixelized_dict={"hello": "there"},
+            search_pixelization_dict={"hello": "there"},
         )
         assert setup.search_inversion_cls == af.DynestyDynamic
-        assert setup.search_pixelized_dict == {"hello": "there"}
+        assert setup.search_pixelization_dict == {"hello": "there"}
 
         setup = ag.SetupHyper(search_noise_cls=None, search_noise_dict=None)
         assert setup.search_noise_cls == af.DynestyStatic
diff --git a/test_autogalaxy/config/general.ini b/test_autogalaxy/config/general.ini
index 89df5ca81..4187d7879 100644
--- a/test_autogalaxy/config/general.ini
+++ b/test_autogalaxy/config/general.ini
@@ -4,6 +4,7 @@ log_file=output.log
 log_level=INFO
 samples_to_csv=False
 model_results_decimal_places=3
+info_whitespace_length=80
 remove_files=False
 force_pickle_overwrite=False
 identifier_version=4
diff --git a/test_autogalaxy/config/lensing/priors/test_analysis.json b/test_autogalaxy/config/lensing/priors/test_analysis.json
index 16bac7a1e..32fac4991 100644
--- a/test_autogalaxy/config/lensing/priors/test_analysis.json
+++ b/test_autogalaxy/config/lensing/priors/test_analysis.json
@@ -1,5 +1,5 @@
 {
-    "aa.m.MockPixelization": {
+    "aa.m.MockMesh": {
         "number_clusters": {
             "type": "Uniform",
             "lower_limit": 1.0,
diff --git a/test_autogalaxy/config/priors/pixelizations.json b/test_autogalaxy/config/priors/mesh.json
similarity index 100%
rename from test_autogalaxy/config/priors/pixelizations.json
rename to test_autogalaxy/config/priors/mesh.json
diff --git a/test_autogalaxy/config/visualize/include.ini b/test_autogalaxy/config/visualize/include.ini
index 9b212dce2..dfbeec7db 100644
--- a/test_autogalaxy/config/visualize/include.ini
+++ b/test_autogalaxy/config/visualize/include.ini
@@ -12,8 +12,8 @@ caustics=False
 multiple_images=False
 light_profile_centres=True
 mass_profile_centres=True
-mapper_data_pixelization_grid=True
-mapper_source_pixelization_grid=True
+mapper_data_mesh_grid=True
+mapper_source_mesh_grid=False
 mapper_source_grid_slim=False
 parallel_overscan=True
 serial_prescan=True
diff --git a/test_autogalaxy/galaxy/test_dict.py b/test_autogalaxy/galaxy/test_dict.py
index 4f7654c35..9c7f4ef2d 100644
--- a/test_autogalaxy/galaxy/test_dict.py
+++ b/test_autogalaxy/galaxy/test_dict.py
@@ -1,7 +1,7 @@
 import pytest
 
 import autogalaxy as ag
-from autoarray.inversion.pixelizations.voronoi import Voronoi
+from autoarray.inversion.pixelization.mesh.voronoi import Voronoi
 from autoarray.inversion.regularization import AdaptiveBrightness
 from autoconf.dictable import Dictable
 
@@ -15,30 +15,32 @@ def make_trivial_galaxy():
 def make_trivial_galaxy_dict():
     return {
         "hyper_galaxy": None,
-        "pixelization": None,
         "redshift": 1.0,
-        "regularization": None,
         "type": "autogalaxy.galaxy.galaxy.Galaxy",
     }
 
 
-class TestTrivial:
-    def test_dict(self, trivial_galaxy, trivial_galaxy_dict):
-        assert trivial_galaxy.dict() == trivial_galaxy_dict
+def test__trivial_dict(trivial_galaxy, trivial_galaxy_dict):
 
-    def test_from_dict(self, trivial_galaxy, trivial_galaxy_dict):
-        assert Dictable.from_dict(trivial_galaxy_dict) == trivial_galaxy
+    assert trivial_galaxy.dict() == trivial_galaxy_dict
 
 
-@pytest.fixture(name="ization_galaxy")
-def make_ization_galaxy():
+def test_trivial__from_dict(trivial_galaxy, trivial_galaxy_dict):
+    assert Dictable.from_dict(trivial_galaxy_dict) == trivial_galaxy
+
+
+@pytest.fixture(name="pixelization_galaxy")
+def make_pixelization_galaxy():
     return ag.Galaxy(
-        redshift=1.0, pixelization=Voronoi(), regularization=AdaptiveBrightness()
+        redshift=1.0,
+        pixelization=ag.Pixelization(
+            mesh=Voronoi(), regularization=AdaptiveBrightness()
+        ),
     )
 
 
-@pytest.fixture(name="ization_galaxy_dict")
-def make_ization_galaxy_dict():
+@pytest.fixture(name="pixelization_galaxy_dict")
+def make_pixelization_galaxy_dict():
     return {
         "hyper_galaxy": None,
         "pixelization": {"type": "autoarray.inversion.pixelizations.voronoi.Voronoi"},
@@ -53,13 +55,17 @@ def make_ization_galaxy_dict():
     }
 
 
-class TestIzations:
-    def test_dict(self, ization_galaxy, ization_galaxy_dict):
-        assert ization_galaxy.dict() == ization_galaxy_dict
+# TODO : Rich fix these to work with Pixelization objects (e.g. contained of mesh and regularization)
+
+
+def test__with_pixelization__dict(pixelization_galaxy, pixelization_galaxy_dict):
 
-    def test_from_dict(self, ization_galaxy, ization_galaxy_dict):
-        galaxy = Dictable.from_dict(ization_galaxy_dict)
-        assert galaxy == ization_galaxy
+    assert pixelization_galaxy.dict() == pixelization_galaxy_dict
+
+
+def test__with_pixelization__from_dict(pixelization_galaxy, pixelization_galaxy_dict):
+    galaxy = Dictable.from_dict(pixelization_galaxy_dict)
+    assert galaxy == pixelization_galaxy
 
 
 def test_pixelization_equality():
@@ -93,16 +99,14 @@ def make_profiles_galaxy_dict():
             "mass_to_light_ratio": 1.0,
             "type": "autogalaxy.profiles.mass_profiles.stellar_mass_profiles.EllDevVaucouleurs",
         },
-        "pixelization": None,
         "redshift": 2.0,
-        "regularization": None,
         "type": "autogalaxy.galaxy.galaxy.Galaxy",
     }
 
 
-class TestProfiles:
-    def test_dict(self, profiles_galaxy, profiles_galaxy_dict):
-        assert profiles_galaxy.dict() == profiles_galaxy_dict
+def test__profiles_galaxy__dict(profiles_galaxy, profiles_galaxy_dict):
+    assert profiles_galaxy.dict() == profiles_galaxy_dict
+
 
-    def test_from_dict(self, profiles_galaxy, profiles_galaxy_dict):
-        assert Dictable.from_dict(profiles_galaxy_dict) == profiles_galaxy
+def test__profiles_galaxy__from_dict(profiles_galaxy, profiles_galaxy_dict):
+    assert Dictable.from_dict(profiles_galaxy_dict) == profiles_galaxy
diff --git a/test_autogalaxy/galaxy/test_galaxy.py b/test_autogalaxy/galaxy/test_galaxy.py
index c6061920c..7bc79005b 100644
--- a/test_autogalaxy/galaxy/test_galaxy.py
+++ b/test_autogalaxy/galaxy/test_galaxy.py
@@ -1,7 +1,8 @@
-import autogalaxy as ag
 import numpy as np
 import pytest
 
+import autogalaxy as ag
+
 from autogalaxy import exc
 
 
@@ -412,16 +413,6 @@ def test__image_2d_from__does_not_include_linear_light_profiles(sub_grid_2d_7x7,
     assert (image == lp_image).all()
 
 
-def test__only_pixelization_raises_error():
-    with pytest.raises(exc.GalaxyException):
-        ag.Galaxy(redshift=0.5, pixelization=object())
-
-
-def test__only_regularization_raises_error():
-    with pytest.raises(exc.GalaxyException):
-        ag.Galaxy(redshift=0.5, regularization=object())
-
-
 def test__light_profile_2d_quantity_from_grid__symmetric_profiles_give_symmetric_results():
     lp_0 = ag.lp.EllSersic(centre=(0.0, 0.0), intensity=1.0)
 
diff --git a/test_autogalaxy/imaging/model/test_visualizer_imaging.py b/test_autogalaxy/imaging/model/test_visualizer_imaging.py
index 31f65e87e..989ff8a09 100644
--- a/test_autogalaxy/imaging/model/test_visualizer_imaging.py
+++ b/test_autogalaxy/imaging/model/test_visualizer_imaging.py
@@ -19,51 +19,41 @@ def push_config(plot_path):
     conf.instance.push(path.join(directory, "config"), output_path=plot_path)
 
 
-class TestVisualizer:
-    def test__visualizes_fit_imaging__uses_configs(
-        self,
-        masked_imaging_7x7,
-        fit_imaging_x2_galaxy_inversion_7x7,
-        include_2d_all,
-        plot_path,
-        plot_patch,
-    ):
-
-        if path.exists(plot_path):
-            shutil.rmtree(plot_path)
-
-        visualizer = VisualizerImaging(visualize_path=plot_path)
-
-        visualizer.visualize_fit_imaging(
-            fit=fit_imaging_x2_galaxy_inversion_7x7, during_analysis=False
-        )
-
-        plot_path = path.join(plot_path, "fit_imaging")
-
-        assert path.join(plot_path, "subplot_fit_imaging.png") in plot_patch.paths
-        assert path.join(plot_path, "image_2d.png") in plot_patch.paths
-        assert path.join(plot_path, "noise_map.png") not in plot_patch.paths
-        assert path.join(plot_path, "signal_to_noise_map.png") not in plot_patch.paths
-        assert path.join(plot_path, "model_image.png") in plot_patch.paths
-        assert path.join(plot_path, "residual_map.png") not in plot_patch.paths
-        assert path.join(plot_path, "normalized_residual_map.png") in plot_patch.paths
-        assert path.join(plot_path, "chi_squared_map.png") in plot_patch.paths
-        assert (
-            path.join(plot_path, "subtracted_image_of_galaxy_0.png") in plot_patch.paths
-        )
-        assert (
-            path.join(plot_path, "subtracted_image_of_galaxy_1.png") in plot_patch.paths
-        )
-
-        assert (
-            path.join(plot_path, "model_image_of_galaxy_0.png") not in plot_patch.paths
-        )
-        assert (
-            path.join(plot_path, "model_image_of_galaxy_1.png") not in plot_patch.paths
-        )
-
-        image = ag.util.array_2d.numpy_array_2d_via_fits_from(
-            file_path=path.join(plot_path, "fits", "image_2d.fits"), hdu=0
-        )
-
-        assert image.shape == (5, 5)
+def test__visualizes_fit_imaging__uses_configs(
+    masked_imaging_7x7,
+    fit_imaging_x2_galaxy_inversion_7x7,
+    include_2d_all,
+    plot_path,
+    plot_patch,
+):
+
+    if path.exists(plot_path):
+        shutil.rmtree(plot_path)
+
+    visualizer = VisualizerImaging(visualize_path=plot_path)
+
+    visualizer.visualize_fit_imaging(
+        fit=fit_imaging_x2_galaxy_inversion_7x7, during_analysis=False
+    )
+
+    plot_path = path.join(plot_path, "fit_imaging")
+
+    assert path.join(plot_path, "subplot_fit_imaging.png") in plot_patch.paths
+    assert path.join(plot_path, "image_2d.png") in plot_patch.paths
+    assert path.join(plot_path, "noise_map.png") not in plot_patch.paths
+    assert path.join(plot_path, "signal_to_noise_map.png") not in plot_patch.paths
+    assert path.join(plot_path, "model_image.png") in plot_patch.paths
+    assert path.join(plot_path, "residual_map.png") not in plot_patch.paths
+    assert path.join(plot_path, "normalized_residual_map.png") in plot_patch.paths
+    assert path.join(plot_path, "chi_squared_map.png") in plot_patch.paths
+    assert path.join(plot_path, "subtracted_image_of_galaxy_0.png") in plot_patch.paths
+    assert path.join(plot_path, "subtracted_image_of_galaxy_1.png") in plot_patch.paths
+
+    assert path.join(plot_path, "model_image_of_galaxy_0.png") not in plot_patch.paths
+    assert path.join(plot_path, "model_image_of_galaxy_1.png") not in plot_patch.paths
+
+    image = ag.util.array_2d.numpy_array_2d_via_fits_from(
+        file_path=path.join(plot_path, "fits", "image_2d.fits"), hdu=0
+    )
+
+    assert image.shape == (5, 5)
diff --git a/test_autogalaxy/imaging/test_fit_imaging.py b/test_autogalaxy/imaging/test_fit_imaging.py
index ded084334..e3662d164 100644
--- a/test_autogalaxy/imaging/test_fit_imaging.py
+++ b/test_autogalaxy/imaging/test_fit_imaging.py
@@ -10,9 +10,7 @@ def test__model_image__with_and_without_psf_blurring(
 
     g0 = ag.Galaxy(
         redshift=0.5,
-        light_profile=ag.m.MockLightProfile(
-            image_2d_value=1.0, image_2d_first_value=2.0
-        ),
+        bulge=ag.m.MockLightProfile(image_2d_value=1.0, image_2d_first_value=2.0),
     )
     plane = ag.Plane(galaxies=[g0])
 
@@ -33,9 +31,7 @@ def test__model_image__with_and_without_psf_blurring(
 
 def test__noise_map__with_and_without_hyper_galaxy(masked_imaging_7x7_no_blur):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d_value=1.0)
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d_value=1.0))
 
     plane = ag.Plane(galaxies=[g0])
 
@@ -49,7 +45,7 @@ def test__noise_map__with_and_without_hyper_galaxy(masked_imaging_7x7_no_blur):
 
     g0 = ag.Galaxy(
         redshift=0.5,
-        light_profile=ag.m.MockLightProfile(image_2d_value=1.0),
+        bulge=ag.m.MockLightProfile(image_2d_value=1.0),
         hyper_galaxy=ag.HyperGalaxy(
             contribution_factor=1.0, noise_factor=1.0, noise_power=1.0
         ),
@@ -74,7 +70,7 @@ def test__noise_map__with_hyper_galaxy_reaches_upper_limit(masked_imaging_7x7_no
 
     g0 = ag.Galaxy(
         redshift=0.5,
-        light_profile=ag.m.MockLightProfile(image_2d_value=1.0),
+        bulge=ag.m.MockLightProfile(image_2d_value=1.0),
         hyper_galaxy=ag.HyperGalaxy(
             contribution_factor=1.0, noise_factor=1.0e9, noise_power=1.0
         ),
@@ -95,9 +91,7 @@ def test__noise_map__with_hyper_galaxy_reaches_upper_limit(masked_imaging_7x7_no
 
 def test__image__with_and_without_hyper_background_sky(masked_imaging_7x7_no_blur):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d_value=1.0)
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d_value=1.0))
 
     plane = ag.Plane(galaxies=[g0])
 
@@ -117,9 +111,7 @@ def test__image__with_and_without_hyper_background_sky(masked_imaging_7x7_no_blu
 
 def test__noise_map__with_and_without_hyper_background(masked_imaging_7x7_no_blur):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d_value=1.0)
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d_value=1.0))
     plane = ag.Plane(galaxies=[g0])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7_no_blur, plane=plane)
@@ -144,71 +136,122 @@ def test__noise_map__with_and_without_hyper_background(masked_imaging_7x7_no_blu
 
 def test__fit_figure_of_merit(masked_imaging_7x7):
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_likelihood == pytest.approx(-75938.05, 1e-4)
+    assert fit.perform_inversion is False
     assert fit.figure_of_merit == pytest.approx(-75938.05, 1.0e-4)
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp.EllSersic(intensity=1.0),
+            ag.lp.EllSersic(intensity=1.0),
+        ]
+    )
 
-    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(redshift=0.5, galaxies=[g0])
+
+    fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
+
+    assert fit.perform_inversion is False
+    assert fit.figure_of_merit == pytest.approx(-75938.05, 1.0e-4)
+
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
+
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), galaxy_pix])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_evidence == pytest.approx(-22.9005, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-22.9005, 1.0e-4)
 
-    galaxy_light = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    galaxy_light = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_pix])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_evidence == pytest.approx(-6840.5851, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-6840.5851, 1.0e-4)
 
     g0_linear_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear.EllSersic(sersic_index=1.0)
+        redshift=0.5, bulge=ag.lp_linear.EllSersic(sersic_index=1.0)
     )
 
     g1_linear_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear.EllSersic(sersic_index=4.0)
+        redshift=0.5, bulge=ag.lp_linear.EllSersic(sersic_index=4.0)
     )
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, g1_linear_light])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_likelihood == pytest.approx(-14.52327, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-14.52327, 1.0e-4)
 
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp_linear.EllSersic(sersic_index=1.0),
+            ag.lp_linear.EllSersic(sersic_index=4.0),
+        ]
+    )
+
+    g0 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(redshift=0.5, galaxies=[g0])
+
+    fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
+
+    assert fit.perform_inversion is True
+    assert fit.figure_of_merit == pytest.approx(-14.52327, 1.0e-4)
+
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp_linear.EllSersic(sersic_index=1.0),
+            ag.lp_linear.EllSersic(sersic_index=4.0),
+        ],
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
+
+    g0 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(redshift=0.5, galaxies=[g0])
+
+    fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
+
+    assert fit.perform_inversion is True
+    assert fit.figure_of_merit == pytest.approx(-29.21448984, 1.0e-4)
+
     g0_operated_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_operated.EllSersic(intensity=1.0)
+        redshift=0.5, bulge=ag.lp_operated.EllSersic(intensity=1.0)
     )
     g1_operated_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_operated.EllSersic(intensity=1.0)
+        redshift=0.5, bulge=ag.lp_operated.EllSersic(intensity=1.0)
     )
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_operated_light, g1_operated_light])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_likelihood == pytest.approx(-342374.9618, 1e-4)
+    assert fit.perform_inversion is False
     assert fit.figure_of_merit == pytest.approx(-342374.9618, 1.0e-4)
 
     g0_linear_operated_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear_operated.EllSersic(sersic_index=1.0)
+        redshift=0.5, bulge=ag.lp_linear_operated.EllSersic(sersic_index=1.0)
     )
     g1_linear_operated_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear_operated.EllSersic(sersic_index=4.0)
+        redshift=0.5, bulge=ag.lp_linear_operated.EllSersic(sersic_index=4.0)
     )
 
     plane = ag.Plane(
@@ -217,14 +260,14 @@ def test__fit_figure_of_merit(masked_imaging_7x7):
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_likelihood == pytest.approx(-14.7237273, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-14.7237273, 1.0e-4)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, galaxy_pix])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.log_evidence == pytest.approx(-22.87827302, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-22.87827302, 1.0e-4)
 
 
@@ -239,13 +282,13 @@ def test__fit_figure_of_merit__include_hyper_methods(masked_imaging_7x7):
 
     g0 = ag.Galaxy(
         redshift=0.5,
-        light_profile=ag.lp.EllSersic(intensity=1.0),
+        bulge=ag.lp.EllSersic(intensity=1.0),
         hyper_galaxy=hyper_galaxy,
         hyper_model_image=np.ones(9),
         hyper_galaxy_image=np.ones(9),
         hyper_minimum_value=0.0,
     )
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
@@ -258,16 +301,16 @@ def test__fit_figure_of_merit__include_hyper_methods(masked_imaging_7x7):
 
     assert (fit.image == np.full(fill_value=2.0, shape=(9,))).all()
     assert (fit.noise_map == np.full(fill_value=5.0, shape=(9,))).all()
-    assert fit.log_likelihood == pytest.approx(-12104.68, 1e-4)
     assert fit.figure_of_merit == pytest.approx(-12104.68, 1.0e-4)
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
     g0 = ag.Galaxy(
         redshift=0.5,
-        pixelization=pix,
-        regularization=reg,
+        pixelization=pixelization,
         hyper_galaxy=hyper_galaxy,
         hyper_model_image=np.ones(9),
         hyper_galaxy_image=np.ones(9),
@@ -286,19 +329,18 @@ def test__fit_figure_of_merit__include_hyper_methods(masked_imaging_7x7):
 
     assert (fit.image == np.full(fill_value=2.0, shape=(9,))).all()
     assert (fit.noise_map == np.full(fill_value=5.0, shape=(9,))).all()
-    assert fit.log_evidence == pytest.approx(-30.1448, 1e-4)
     assert fit.figure_of_merit == pytest.approx(-30.1448, 1.0e-4)
 
     galaxy_light = ag.Galaxy(
         redshift=0.5,
-        light_profile=ag.lp.EllSersic(intensity=1.0),
+        bulge=ag.lp.EllSersic(intensity=1.0),
         hyper_galaxy=hyper_galaxy,
         hyper_model_image=np.ones(9),
         hyper_galaxy_image=np.ones(9),
         hyper_minimum_value=0.0,
     )
 
-    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_pix])
 
@@ -312,7 +354,6 @@ def test__fit_figure_of_merit__include_hyper_methods(masked_imaging_7x7):
 
     assert (fit.image == np.full(fill_value=2.0, shape=(9,))).all()
     assert (fit.noise_map == np.full(fill_value=5.0, shape=(9,))).all()
-    assert fit.log_evidence == pytest.approx(-1132.9297, 1e-4)
     assert fit.figure_of_merit == pytest.approx(-1132.9297, 1.0e-4)
 
 
@@ -320,8 +361,8 @@ def test__galaxy_model_image_dict(masked_imaging_7x7):
 
     # Normal Light Profiles Only
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=2.0))
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=2.0))
     g2 = ag.Galaxy(
         redshift=0.5,
         light_profile_0=ag.lp.EllSersic(intensity=1.0),
@@ -361,7 +402,7 @@ def test__galaxy_model_image_dict(masked_imaging_7x7):
 
     # Linear Light Profiles only
 
-    g0_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    g0_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear, g3])
 
@@ -378,11 +419,13 @@ def test__galaxy_model_image_dict(masked_imaging_7x7):
 
     # Pixelization + Regularizaiton only
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
     g0 = ag.Galaxy(redshift=0.5)
-    g1 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
@@ -401,14 +444,16 @@ def test__galaxy_model_image_dict(masked_imaging_7x7):
 
     # Linear Light PRofiles + Pixelization + Regularizaiton
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    g1_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    g1_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
-    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
-    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
+    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(
         redshift=0.5, galaxies=[g0, g1_linear, g3, galaxy_pix_0, galaxy_pix_1]
@@ -449,15 +494,16 @@ def test__galaxy_model_image_dict(masked_imaging_7x7):
 
 def test__model_images_of_galaxies_list(masked_imaging_7x7):
 
-    galaxy_light = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    galaxy_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    galaxy_light = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    galaxy_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.pix.Rectangular(shape=(3, 3)),
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
         regularization=ag.reg.Constant(coefficient=1.0),
     )
 
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
+
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_linear, galaxy_pix])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
@@ -474,9 +520,9 @@ def test__model_images_of_galaxies_list(masked_imaging_7x7):
 
 
 def test___unmasked_blurred_images(masked_imaging_7x7):
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
@@ -504,17 +550,11 @@ def test___unmasked_blurred_images(masked_imaging_7x7):
 
 def test__subtracted_images_of_galaxies(masked_imaging_7x7_no_blur):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=np.ones(1))
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=np.ones(1)))
 
-    g1 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=2.0 * np.ones(1))
-    )
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=2.0 * np.ones(1)))
 
-    g2 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=3.0 * np.ones(1))
-    )
+    g2 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=3.0 * np.ones(1)))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1, g2])
 
@@ -526,15 +566,11 @@ def test__subtracted_images_of_galaxies(masked_imaging_7x7_no_blur):
     assert fit.subtracted_images_of_galaxies_list[1].slim[0] == -3.0 or np.nan
     assert fit.subtracted_images_of_galaxies_list[2].slim[0] == -2.0 or np.nan
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=np.ones(1))
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=np.ones(1)))
 
     g1 = ag.Galaxy(redshift=0.5)
 
-    g2 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=3.0 * np.ones(1))
-    )
+    g2 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=3.0 * np.ones(1)))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1, g2])
 
@@ -550,9 +586,9 @@ def test__light_profile_linear__intensity_dict(masked_imaging_7x7):
     linear_light_0 = ag.lp_linear.EllSersic(sersic_index=1.0)
     linear_light_1 = ag.lp_linear.EllSersic(sersic_index=4.0)
 
-    g0_linear_light = ag.Galaxy(redshift=0.5, light_profile=linear_light_0)
+    g0_linear_light = ag.Galaxy(redshift=0.5, bulge=linear_light_0)
 
-    g1_linear_light = ag.Galaxy(redshift=0.5, light_profile=linear_light_1)
+    g1_linear_light = ag.Galaxy(redshift=0.5, bulge=linear_light_1)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, g1_linear_light])
 
@@ -565,42 +601,72 @@ def test__light_profile_linear__intensity_dict(masked_imaging_7x7):
         -0.04694839915145, 1.0e-4
     )
 
+    basis = ag.lp_basis.Basis(light_profile_list=[linear_light_0, linear_light_1])
 
-def test__plane_linear_light_profiles_to_light_profiles(masked_imaging_7x7):
+    g_basis = ag.Galaxy(redshift=0.5, bulge=basis)
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    plane = ag.Plane(redshift=0.5, galaxies=[g_basis])
 
-    g0_linear = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear.EllSersic(sersic_index=1.0)
+    fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
+
+    assert fit.linear_light_profile_intensity_dict[linear_light_0] == pytest.approx(
+        7.093227476666252, 1.0e-4
     )
+    assert fit.linear_light_profile_intensity_dict[linear_light_1] == pytest.approx(
+        -0.04694839915145, 1.0e-4
+    )
+
+    linear_light_2 = ag.lp_linear.EllSersic(sersic_index=2.0)
+    linear_light_3 = ag.lp_linear.EllSersic(sersic_index=3.0)
+
+    basis = ag.lp_basis.Basis(light_profile_list=[linear_light_2, linear_light_3])
+
+    g_basis = ag.Galaxy(redshift=0.5, bulge=basis)
 
-    g1_linear = ag.Galaxy(
-        redshift=1.0, light_profile=ag.lp_linear.EllSersic(sersic_index=4.0)
+    plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, g_basis])
+
+    fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
+
+    assert fit.linear_light_profile_intensity_dict[linear_light_0] == pytest.approx(
+        -21.77759470, 1.0e-4
     )
+    assert fit.linear_light_profile_intensity_dict[linear_light_2] == pytest.approx(
+        29.3935231947, 1.0e-4
+    )
+    assert fit.linear_light_profile_intensity_dict[linear_light_3] == pytest.approx(
+        -4.77469646, 1.0e-4
+    )
+
+
+def test__plane_linear_light_profiles_to_light_profiles(masked_imaging_7x7):
+
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+
+    g0_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic(sersic_index=1.0))
+
+    g1_linear = ag.Galaxy(redshift=1.0, bulge=ag.lp_linear.EllSersic(sersic_index=4.0))
 
     plane = ag.Plane(galaxies=[g0, g0_linear, g1_linear])
 
     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 
-    assert fit.galaxies[0].light_profile.intensity == pytest.approx(1.0, 1.0e-4)
+    assert fit.galaxies[0].bulge.intensity == pytest.approx(1.0, 1.0e-4)
 
     plane = fit.plane_linear_light_profiles_to_light_profiles
 
-    assert plane.galaxies[0].light_profile.intensity == pytest.approx(1.0, 1.0e-4)
-    assert plane.galaxies[1].light_profile.intensity == pytest.approx(7.0932274, 1.0e-4)
-    assert plane.galaxies[2].light_profile.intensity == pytest.approx(
-        -1.04694839, 1.0e-4
-    )
+    assert plane.galaxies[0].bulge.intensity == pytest.approx(1.0, 1.0e-4)
+    assert plane.galaxies[1].bulge.intensity == pytest.approx(7.0932274, 1.0e-4)
+    assert plane.galaxies[2].bulge.intensity == pytest.approx(-1.04694839, 1.0e-4)
 
 
 # def test__light_profile_no_convolve(masked_imaging_7x7):
 #
-#     g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllGaussian(intensity=1.0))
+#     g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllGaussian(intensity=1.0))
 #     plane = ag.Plane(redshift=0.5, galaxies=[g0])
 #     fit = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane)
 #
 #     g0_no_convolve = ag.Galaxy(
-#         redshift=0.5, light_profile=ag.lp.EllGaussian(intensity=1.0)
+#         redshift=0.5, bulge=ag.lp.EllGaussian(intensity=1.0)
 #     )
 #     plane_no_convolve = ag.Plane(redshift=0.5, galaxies=[g0_no_convolve])
 #     fit_no_convolve = ag.FitImaging(dataset=masked_imaging_7x7, plane=plane_no_convolve)
diff --git a/test_autogalaxy/imaging/test_imaging.py b/test_autogalaxy/imaging/test_imaging.py
index 4395f4087..68b4df373 100644
--- a/test_autogalaxy/imaging/test_imaging.py
+++ b/test_autogalaxy/imaging/test_imaging.py
@@ -31,128 +31,119 @@ def clean_fits(fits_path):
         shutil.rmtree(fits_path)
 
 
-class TestImaging:
-    def test__from_fits__all_imaging_data_structures_are_flipped_for_ds9(self):
+def test__from_fits__all_imaging_data_structures_are_flipped_for_ds9():
 
-        fits_path = path.join(
-            "{}".format(path.dirname(path.realpath(__file__))), "files"
-        )
-        image_path = path.join(fits_path, "image.fits")
-        noise_map_path = path.join(fits_path, "noise_map.fits")
-        psf_path = path.join(fits_path, "psf.fits")
-
-        create_fits(fits_path=image_path, array=[[1.0, 0.0], [0.0, 0.0]])
-        create_fits(fits_path=noise_map_path, array=[[2.0, 0.0], [0.0, 0.0]])
-        create_fits(fits_path=psf_path, array=[[1.0, 1.0], [0.0, 0.0]])
-
-        imaging = ag.Imaging.from_fits(
-            image_path=image_path,
-            noise_map_path=noise_map_path,
-            psf_path=psf_path,
-            pixel_scales=0.1,
-        )
+    fits_path = path.join("{}".format(path.dirname(path.realpath(__file__))), "files")
+    image_path = path.join(fits_path, "image.fits")
+    noise_map_path = path.join(fits_path, "noise_map.fits")
+    psf_path = path.join(fits_path, "psf.fits")
 
-        assert (imaging.image.native == np.array([[0.0, 0.0], [1.0, 0.0]])).all()
-        assert (imaging.noise_map.native == np.array([[0.0, 0.0], [2.0, 0.0]])).all()
-        assert (imaging.psf.native == np.array([[0.0, 0.0], [0.5, 0.5]])).all()
-
-        imaging.output_to_fits(
-            image_path=image_path,
-            noise_map_path=noise_map_path,
-            psf_path=psf_path,
-            overwrite=True,
-        )
-
-        hdu_list = fits.open(image_path)
-        image = np.array(hdu_list[0].data).astype("float64")
-        assert (image == np.array([[1.0, 0.0], [0.0, 0.0]])).all()
-
-        hdu_list = fits.open(noise_map_path)
-        noise_map = np.array(hdu_list[0].data).astype("float64")
-        assert (noise_map == np.array([[2.0, 0.0], [0.0, 0.0]])).all()
-
-        hdu_list = fits.open(psf_path)
-        psf = np.array(hdu_list[0].data).astype("float64")
-        assert (psf == np.array([[0.5, 0.5], [0.0, 0.0]])).all()
-
-        clean_fits(fits_path=fits_path)
-
-
-class TestSimulatorImaging:
-    def test__via_plane_from__same_as_plane_image(self):
-
-        psf = ag.Kernel2D.from_gaussian(
-            shape_native=(7, 7), sigma=0.5, pixel_scales=0.05
-        )
-
-        grid = ag.Grid2D.uniform(shape_native=(20, 20), pixel_scales=0.05, sub_size=1)
-
-        galaxy_0 = ag.Galaxy(redshift=0.5, light=ag.lp.EllSersic(intensity=1.0))
-
-        galaxy_1 = ag.Galaxy(redshift=1.0, light=ag.lp.EllSersic(intensity=0.3))
-
-        plane = ag.Plane(redshift=0.75, galaxies=[galaxy_0, galaxy_1])
-
-        simulator = ag.SimulatorImaging(
-            psf=psf,
-            exposure_time=10000.0,
-            background_sky_level=100.0,
-            add_poisson_noise=False,
-        )
-
-        imaging = simulator.via_plane_from(plane=plane, grid=grid)
-
-        imaging_via_image = simulator.via_image_from(
-            image=plane.image_2d_from(grid=grid)
-        )
-
-        assert imaging.shape_native == (20, 20)
-        assert imaging.image.native[0, 0] != imaging_via_image.image.native[0, 0]
-        assert imaging.image.native[10, 10] == imaging_via_image.image.native[10, 10]
-        assert (imaging.psf == imaging_via_image.psf).all()
-        assert (imaging.noise_map == imaging_via_image.noise_map).all()
-
-    def test__simulate_imaging_from_galaxy__source_galaxy__compare_to_imaging(self):
-
-        galaxy_0 = ag.Galaxy(
-            redshift=0.5,
-            mass=ag.mp.EllIsothermal(
-                centre=(0.0, 0.0), einstein_radius=1.6, elliptical_comps=(0.17647, 0.0)
-            ),
-        )
-
-        galaxy_1 = ag.Galaxy(
-            redshift=1.0,
-            light=ag.lp.EllSersic(
-                centre=(0.1, 0.1),
-                elliptical_comps=(0.096225, -0.055555),
-                intensity=0.3,
-                effective_radius=1.0,
-                sersic_index=2.5,
-            ),
-        )
-
-        grid = ag.Grid2D.uniform(shape_native=(11, 11), pixel_scales=0.2, sub_size=1)
+    create_fits(fits_path=image_path, array=[[1.0, 0.0], [0.0, 0.0]])
+    create_fits(fits_path=noise_map_path, array=[[2.0, 0.0], [0.0, 0.0]])
+    create_fits(fits_path=psf_path, array=[[1.0, 1.0], [0.0, 0.0]])
 
-        psf = ag.Kernel2D.manual_native(array=[[1.0]], pixel_scales=0.2)
-
-        simulator = ag.SimulatorImaging(
-            psf=psf,
-            exposure_time=10000.0,
-            background_sky_level=100.0,
-            add_poisson_noise=True,
-            noise_seed=1,
-        )
+    imaging = ag.Imaging.from_fits(
+        image_path=image_path,
+        noise_map_path=noise_map_path,
+        psf_path=psf_path,
+        pixel_scales=0.1,
+    )
 
-        imaging = simulator.via_galaxies_from(galaxies=[galaxy_0, galaxy_1], grid=grid)
+    assert (imaging.image.native == np.array([[0.0, 0.0], [1.0, 0.0]])).all()
+    assert (imaging.noise_map.native == np.array([[0.0, 0.0], [2.0, 0.0]])).all()
+    assert (imaging.psf.native == np.array([[0.0, 0.0], [0.5, 0.5]])).all()
 
-        plane = ag.Plane(redshift=0.75, galaxies=[galaxy_0, galaxy_1])
+    imaging.output_to_fits(
+        image_path=image_path,
+        noise_map_path=noise_map_path,
+        psf_path=psf_path,
+        overwrite=True,
+    )
 
-        imaging_via_image = simulator.via_image_from(
-            image=plane.image_2d_from(grid=grid)
-        )
+    hdu_list = fits.open(image_path)
+    image = np.array(hdu_list[0].data).astype("float64")
+    assert (image == np.array([[1.0, 0.0], [0.0, 0.0]])).all()
 
-        assert imaging.shape_native == (11, 11)
-        assert imaging.image == pytest.approx(imaging_via_image.image, 1.0e-4)
-        assert (imaging.psf == imaging_via_image.psf).all()
-        assert (imaging.noise_map == imaging_via_image.noise_map).all()
+    hdu_list = fits.open(noise_map_path)
+    noise_map = np.array(hdu_list[0].data).astype("float64")
+    assert (noise_map == np.array([[2.0, 0.0], [0.0, 0.0]])).all()
+
+    hdu_list = fits.open(psf_path)
+    psf = np.array(hdu_list[0].data).astype("float64")
+    assert (psf == np.array([[0.5, 0.5], [0.0, 0.0]])).all()
+
+    clean_fits(fits_path=fits_path)
+
+
+def test__simulator__via_plane_from__same_as_plane_image():
+
+    psf = ag.Kernel2D.from_gaussian(shape_native=(7, 7), sigma=0.5, pixel_scales=0.05)
+
+    grid = ag.Grid2D.uniform(shape_native=(20, 20), pixel_scales=0.05, sub_size=1)
+
+    galaxy_0 = ag.Galaxy(redshift=0.5, light=ag.lp.EllSersic(intensity=1.0))
+
+    galaxy_1 = ag.Galaxy(redshift=1.0, light=ag.lp.EllSersic(intensity=0.3))
+
+    plane = ag.Plane(redshift=0.75, galaxies=[galaxy_0, galaxy_1])
+
+    simulator = ag.SimulatorImaging(
+        psf=psf,
+        exposure_time=10000.0,
+        background_sky_level=100.0,
+        add_poisson_noise=False,
+    )
+
+    imaging = simulator.via_plane_from(plane=plane, grid=grid)
+
+    imaging_via_image = simulator.via_image_from(image=plane.image_2d_from(grid=grid))
+
+    assert imaging.shape_native == (20, 20)
+    assert imaging.image.native[0, 0] != imaging_via_image.image.native[0, 0]
+    assert imaging.image.native[10, 10] == imaging_via_image.image.native[10, 10]
+    assert (imaging.psf == imaging_via_image.psf).all()
+    assert (imaging.noise_map == imaging_via_image.noise_map).all()
+
+
+def test__simulator__simulate_imaging_from_galaxy__source_galaxy__compare_to_imaging():
+
+    galaxy_0 = ag.Galaxy(
+        redshift=0.5,
+        mass=ag.mp.EllIsothermal(
+            centre=(0.0, 0.0), einstein_radius=1.6, elliptical_comps=(0.17647, 0.0)
+        ),
+    )
+
+    galaxy_1 = ag.Galaxy(
+        redshift=1.0,
+        light=ag.lp.EllSersic(
+            centre=(0.1, 0.1),
+            elliptical_comps=(0.096225, -0.055555),
+            intensity=0.3,
+            effective_radius=1.0,
+            sersic_index=2.5,
+        ),
+    )
+
+    grid = ag.Grid2D.uniform(shape_native=(11, 11), pixel_scales=0.2, sub_size=1)
+
+    psf = ag.Kernel2D.manual_native(array=[[1.0]], pixel_scales=0.2)
+
+    simulator = ag.SimulatorImaging(
+        psf=psf,
+        exposure_time=10000.0,
+        background_sky_level=100.0,
+        add_poisson_noise=True,
+        noise_seed=1,
+    )
+
+    imaging = simulator.via_galaxies_from(galaxies=[galaxy_0, galaxy_1], grid=grid)
+
+    plane = ag.Plane(redshift=0.75, galaxies=[galaxy_0, galaxy_1])
+
+    imaging_via_image = simulator.via_image_from(image=plane.image_2d_from(grid=grid))
+
+    assert imaging.shape_native == (11, 11)
+    assert imaging.image == pytest.approx(imaging_via_image.image, 1.0e-4)
+    assert (imaging.psf == imaging_via_image.psf).all()
+    assert (imaging.noise_map == imaging_via_image.noise_map).all()
diff --git a/test_autogalaxy/imaging/test_simulate_and_fit_imaging.py b/test_autogalaxy/imaging/test_simulate_and_fit_imaging.py
index 4afe24cd0..2302832f7 100644
--- a/test_autogalaxy/imaging/test_simulate_and_fit_imaging.py
+++ b/test_autogalaxy/imaging/test_simulate_and_fit_imaging.py
@@ -94,19 +94,23 @@ def test__simulate_imaging_data_and_fit__known_likelihood():
     )
 
     galaxy_0 = ag.Galaxy(
-        redshift=0.5, light=ag.lp.EllSersic(centre=(0.1, 0.1), intensity=0.1)
+        redshift=0.5,
+        bulge=ag.lp.EllSersic(centre=(0.1, 0.1), intensity=0.1),
+        disk=ag.lp.EllSersic(centre=(0.2, 0.2), intensity=0.2),
     )
-    galaxy_1 = ag.Galaxy(
-        redshift=1.0,
-        pixelization=ag.pix.Rectangular(shape=(16, 16)),
+
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(16, 16)),
         regularization=ag.reg.Constant(coefficient=(1.0)),
     )
 
+    galaxy_1 = ag.Galaxy(redshift=1.0, pixelization=pixelization)
+
     plane = ag.Plane(galaxies=[galaxy_0, galaxy_1])
 
-    imaging = ag.SimulatorImaging(exposure_time=300.0, psf=psf, noise_seed=1)
+    simulator = ag.SimulatorImaging(exposure_time=300.0, psf=psf, noise_seed=1)
 
-    imaging = imaging.via_plane_from(plane=plane, grid=grid)
+    imaging = simulator.via_plane_from(plane=plane, grid=grid)
 
     mask = ag.Mask2D.circular(
         shape_native=imaging.image.shape_native, pixel_scales=0.2, radius=2.0
@@ -116,7 +120,28 @@ def test__simulate_imaging_data_and_fit__known_likelihood():
 
     fit = ag.FitImaging(dataset=masked_imaging, plane=plane)
 
-    assert fit.figure_of_merit == pytest.approx(621.5594864161, 1.0e-2)
+    assert fit.figure_of_merit == pytest.approx(532.19918562, 1.0e-2)
+
+    # Check that using a Basis gives the same result.
+
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp.EllSersic(centre=(0.1, 0.1), intensity=0.1),
+            ag.lp.EllSersic(centre=(0.2, 0.2), intensity=0.2),
+        ]
+    )
+
+    galaxy_0 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(galaxies=[galaxy_0, galaxy_1])
+
+    imaging = simulator.via_plane_from(plane=plane, grid=grid)
+
+    masked_imaging = imaging.apply_mask(mask=mask)
+
+    fit = ag.FitImaging(dataset=masked_imaging, plane=plane)
+
+    assert fit.figure_of_merit == pytest.approx(532.19918562, 1.0e-2)
 
 
 def test__simulate_imaging_data_and_fit__linear_light_profiles_agree_with_standard_light_profiles():
diff --git a/test_autogalaxy/interferometer/test_fit_interferometer.py b/test_autogalaxy/interferometer/test_fit_interferometer.py
index 63a879ce5..95ed4e059 100644
--- a/test_autogalaxy/interferometer/test_fit_interferometer.py
+++ b/test_autogalaxy/interferometer/test_fit_interferometer.py
@@ -6,9 +6,7 @@
 
 def test__model_visibilities(interferometer_7):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=np.ones(2))
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=np.ones(2)))
     plane = ag.Plane(galaxies=[g0])
 
     fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
@@ -21,9 +19,7 @@ def test__model_visibilities(interferometer_7):
 
 def test__noise_map__with_and_without_hyper_background(interferometer_7):
 
-    g0 = ag.Galaxy(
-        redshift=0.5, light_profile=ag.m.MockLightProfile(image_2d=np.ones(2))
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.m.MockLightProfile(image_2d=np.ones(2)))
     plane = ag.Plane(galaxies=[g0])
 
     fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
@@ -44,26 +40,39 @@ def test__noise_map__with_and_without_hyper_background(interferometer_7):
 
 def test__fit_figure_of_merit(interferometer_7):
 
-    g0 = ag.Galaxy(
-        redshift=0.5,
-        light_profile=ag.lp.EllSersic(intensity=1.0),
-        mass_profile=ag.mp.SphIsothermal(einstein_radius=1.0),
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
     fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
 
-    assert (fit.noise_map.slim == np.full(fill_value=2.0 + 2.0j, shape=(7,))).all()
-    assert fit.log_likelihood == pytest.approx(-2398107.3849, 1e-4)
+    assert fit.perform_inversion is False
+    assert fit.figure_of_merit == pytest.approx(-2398107.3849, 1.0e-4)
+
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp.EllSersic(intensity=1.0),
+            ag.lp.EllSersic(intensity=1.0),
+        ]
+    )
+
+    g0 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(redshift=0.5, galaxies=[g0])
+
+    fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
+
+    assert fit.perform_inversion is False
     assert fit.figure_of_merit == pytest.approx(-2398107.3849, 1.0e-4)
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=0.01)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=0.01),
+    )
 
-    g0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g0])
 
@@ -73,15 +82,17 @@ def test__fit_figure_of_merit(interferometer_7):
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
-    assert (fit.noise_map.slim == np.full(fill_value=2.0 + 2.0j, shape=(7,))).all()
-    assert fit.log_evidence == pytest.approx(-66.90612, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-66.90612, 1.0e-4)
 
-    galaxy_light = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    galaxy_light = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
-    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_pix])
 
@@ -91,23 +102,36 @@ def test__fit_figure_of_merit(interferometer_7):
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
-    assert (fit.noise_map.slim == np.full(fill_value=2.0 + 2.0j, shape=(7,))).all()
-    assert fit.log_evidence == pytest.approx(-283424.48941, 1e-4)
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-283424.48941, 1.0e-4)
 
     g0_linear_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear.EllSersic(sersic_index=1.0)
+        redshift=0.5, bulge=ag.lp_linear.EllSersic(sersic_index=1.0)
     )
 
     g1_linear_light = ag.Galaxy(
-        redshift=0.5, light_profile=ag.lp_linear.EllSersic(sersic_index=4.0)
+        redshift=0.5, bulge=ag.lp_linear.EllSersic(sersic_index=4.0)
     )
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, g1_linear_light])
 
     fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
 
-    assert fit.log_likelihood == pytest.approx(-23.44419, 1e-4)
+    assert fit.perform_inversion is True
+    assert fit.figure_of_merit == pytest.approx(-23.44419, 1.0e-4)
+
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[
+            ag.lp_linear.EllSersic(sersic_index=1.0),
+            ag.lp_linear.EllSersic(sersic_index=4.0),
+        ]
+    )
+
+    plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, g1_linear_light])
+
+    fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
+
+    assert fit.perform_inversion is True
     assert fit.figure_of_merit == pytest.approx(-23.44419, 1.0e-4)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear_light, galaxy_pix])
@@ -122,13 +146,9 @@ def test__fit_figure_of_merit__include_hyper_methods(interferometer_7):
 
     hyper_background_noise = ag.hyper_data.HyperBackgroundNoise(noise_scale=1.0)
 
-    g0 = ag.Galaxy(
-        redshift=0.5,
-        light_profile=ag.lp.EllSersic(intensity=1.0),
-        mass_profile=ag.mp.SphIsothermal(einstein_radius=1.0),
-    )
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, g1])
 
@@ -153,10 +173,12 @@ def test__fit_figure_of_merit__include_hyper_methods(interferometer_7):
 
     assert fit.noise_map == pytest.approx(interferometer_7.noise_map, 1.0e-4)
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=0.01)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=0.01),
+    )
 
-    g0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g0])
 
@@ -171,11 +193,14 @@ def test__fit_figure_of_merit__include_hyper_methods(interferometer_7):
     assert fit.log_evidence == pytest.approx(-68.63380, 1e-4)
     assert fit.figure_of_merit == pytest.approx(-68.63380, 1.0e-4)
 
-    galaxy_light = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
+    galaxy_light = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
-    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
+
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_pix])
 
@@ -195,10 +220,12 @@ def test___fit_figure_of_merit__different_settings(
     interferometer_7, interferometer_7_lop
 ):
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=0.01)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=0.01),
+    )
 
-    g0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g0])
 
@@ -219,8 +246,8 @@ def test___galaxy_model_image_dict(interferometer_7):
 
     # Normal Light Profiles Only
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=2.0))
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=2.0))
     g2 = ag.Galaxy(
         redshift=0.5,
         light_profile_0=ag.lp.EllSersic(intensity=1.0),
@@ -244,7 +271,7 @@ def test___galaxy_model_image_dict(interferometer_7):
 
     # Linear Light Profiles Only
 
-    g0_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    g0_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear])
 
@@ -260,11 +287,13 @@ def test___galaxy_model_image_dict(interferometer_7):
 
     # Pixelization + Regularizaiton only
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
     g0 = ag.Galaxy(redshift=0.5)
-    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, galaxy_pix_0])
 
@@ -274,14 +303,17 @@ def test___galaxy_model_image_dict(interferometer_7):
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
-    mapper = pix.mapper_from(
-        source_grid_slim=interferometer_7.grid, source_pixelization_grid=None
+    mapper_grids = pixelization.mesh.mapper_grids_from(
+        source_grid_slim=interferometer_7.grid, source_mesh_grid=None
+    )
+
+    mapper = ag.Mapper(
+        mapper_grids=mapper_grids, regularization=pixelization.regularization
     )
 
     inversion = ag.Inversion(
         dataset=interferometer_7,
         linear_obj_list=[mapper],
-        regularization_list=[reg],
         settings=ag.SettingsInversion(use_w_tilde=False),
     )
 
@@ -292,9 +324,13 @@ def test___galaxy_model_image_dict(interferometer_7):
     )
 
     # Linear Light PRofiles + Pixelization + Regularizaiton
-    reg = ag.reg.Constant(coefficient=2.0)
 
-    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=2.0),
+    )
+
+    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear, g1, galaxy_pix_0, galaxy_pix_1])
 
@@ -330,8 +366,8 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
 
     # Normal Light Profiles Only
 
-    g0 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    g1 = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=2.0))
+    g0 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    g1 = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=2.0))
     g2 = ag.Galaxy(
         redshift=0.5,
         light_profile_0=ag.lp.EllSersic(intensity=1.0),
@@ -365,7 +401,7 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
 
     # Linear Light Profiles Only
 
-    g0_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    g0_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear])
 
@@ -381,11 +417,13 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
 
     # Pixelization + Regularizaiton only
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=1.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=1.0),
+    )
 
     g0 = ag.Galaxy(redshift=0.5)
-    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0, galaxy_pix_0])
 
@@ -395,14 +433,17 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
-    mapper = pix.mapper_from(
-        source_grid_slim=interferometer_7.grid, source_pixelization_grid=None
+    mapper_grids = pixelization.mesh.mapper_grids_from(
+        source_grid_slim=interferometer_7.grid, source_mesh_grid=None
+    )
+
+    mapper = ag.Mapper(
+        mapper_grids=mapper_grids, regularization=pixelization.regularization
     )
 
     inversion = ag.Inversion(
         dataset=interferometer_7,
         linear_obj_list=[mapper],
-        regularization_list=[reg],
         settings=ag.SettingsInversion(use_w_tilde=False),
     )
 
@@ -413,9 +454,12 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
     )
 
     # Linear Light PRofiles + Pixelization + Regularizaiton
-    reg = ag.reg.Constant(coefficient=2.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=2.0),
+    )
 
-    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(redshift=0.5, galaxies=[g0_linear, g1, galaxy_pix_0, galaxy_pix_1])
 
@@ -452,15 +496,16 @@ def test___galaxy_model_visibilities_dict(interferometer_7):
 
 def test__model_visibilities_of_galaxies_list(interferometer_7):
 
-    galaxy_light = ag.Galaxy(redshift=0.5, light_profile=ag.lp.EllSersic(intensity=1.0))
-    galaxy_linear = ag.Galaxy(redshift=0.5, light_profile=ag.lp_linear.EllSersic())
+    galaxy_light = ag.Galaxy(redshift=0.5, bulge=ag.lp.EllSersic(intensity=1.0))
+    galaxy_linear = ag.Galaxy(redshift=0.5, bulge=ag.lp_linear.EllSersic())
 
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.pix.Rectangular(shape=(3, 3)),
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
         regularization=ag.reg.Constant(coefficient=1.0),
     )
 
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
+
     plane = ag.Plane(redshift=0.5, galaxies=[galaxy_light, galaxy_linear, galaxy_pix])
 
     fit = ag.FitInterferometer(dataset=interferometer_7, plane=plane)
diff --git a/test_autogalaxy/interferometer/test_simulate_and_fit_interferometer.py b/test_autogalaxy/interferometer/test_simulate_and_fit_interferometer.py
index 4dd5f106b..fa229232d 100644
--- a/test_autogalaxy/interferometer/test_simulate_and_fit_interferometer.py
+++ b/test_autogalaxy/interferometer/test_simulate_and_fit_interferometer.py
@@ -76,15 +76,16 @@ def test__perfect_fit__chi_squared_0():
 
     assert fit.chi_squared == pytest.approx(0.0)
 
-    pix = ag.pix.Rectangular(shape=(7, 7))
-
-    reg = ag.reg.Constant(coefficient=0.0001)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(7, 7)),
+        regularization=ag.reg.Constant(coefficient=0.0001),
+    )
 
     galaxy_0 = ag.Galaxy(
         redshift=0.5, light=ag.lp.EllSersic(centre=(0.1, 0.1), intensity=0.1)
     )
 
-    galaxy_1 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    galaxy_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[galaxy_0, galaxy_1])
 
@@ -116,12 +117,14 @@ def test__simulate_interferometer_data_and_fit__known_likelihood():
         light=ag.lp.EllSersic(centre=(0.1, 0.1), intensity=0.1),
         mass=ag.mp.EllIsothermal(centre=(0.1, 0.1), einstein_radius=1.8),
     )
-    galaxy_1 = ag.Galaxy(
-        redshift=1.0,
-        pixelization=ag.pix.Rectangular(shape=(16, 16)),
+
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(16, 16)),
         regularization=ag.reg.Constant(coefficient=(1.0)),
     )
 
+    galaxy_1 = ag.Galaxy(redshift=1.0, pixelization=pixelization)
+
     plane = ag.Plane(galaxies=[galaxy_0, galaxy_1])
 
     simulator = ag.SimulatorInterferometer(
@@ -199,7 +202,7 @@ def test__linear_light_profiles_agree_with_standard_light_profiles():
         plane=plane_linear,
         settings_pixelization=ag.SettingsPixelization(use_border=False),
         settings_inversion=ag.SettingsInversion(
-            use_w_tilde=False, linear_funcs_add_to_curvature_diag=False
+            use_w_tilde=False, no_regularization_add_to_curvature_diag=False
         ),
     )
 
diff --git a/test_autogalaxy/plane/test_plane.py b/test_autogalaxy/plane/test_plane.py
index 9c2a831bc..385ec185c 100644
--- a/test_autogalaxy/plane/test_plane.py
+++ b/test_autogalaxy/plane/test_plane.py
@@ -4,7 +4,7 @@
 
 import autogalaxy as ag
 
-from autoarray.inversion.pixelizations.abstract import AbstractPixelization
+from autoarray.inversion.pixelization.mesh.abstract import AbstractMesh
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 
 from autogalaxy import exc
@@ -109,54 +109,43 @@ def test__cls_list_from():
     )
     assert plane.cls_list_from(cls=ag.mp.MassProfile) == [sis_0, sis_1, sis_2, sis_1]
 
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(mapper=1),
-        regularization=ag.m.MockRegularization(),
-    )
+    pixelization = ag.m.MockPixelization(mapper=1)
+
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[galaxy_pix], redshift=None)
 
-    assert plane.cls_list_from(cls=ag.pix.Pixelization)[0].mapper == 1
+    assert plane.cls_list_from(cls=ag.Pixelization)[0].mapper == 1
 
-    galaxy_pix_0 = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(mapper=1),
-        regularization=ag.m.MockRegularization(),
-    )
-    galaxy_pix_1 = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(mapper=2),
-        regularization=ag.m.MockRegularization(),
-    )
+    galaxy_pix_0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
+
+    pixelization = ag.m.MockPixelization(mapper=2)
+
+    galaxy_pix_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[galaxy_pix_0, galaxy_pix_1], redshift=None)
 
-    assert plane.cls_list_from(cls=ag.pix.Pixelization)[0].mapper == 1
-    assert plane.cls_list_from(cls=ag.pix.Pixelization)[1].mapper == 2
+    assert plane.cls_list_from(cls=ag.Pixelization)[0].mapper == 1
+    assert plane.cls_list_from(cls=ag.Pixelization)[1].mapper == 2
 
     galaxy_no_pix = ag.Galaxy(redshift=0.5)
 
     plane = ag.Plane(galaxies=[galaxy_no_pix], redshift=None)
 
-    assert plane.cls_list_from(cls=ag.pix.Pixelization) == []
+    assert plane.cls_list_from(cls=ag.Pixelization) == []
 
 
 def test__hyper_galaxy_image_list():
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=AbstractPixelization(),
-        regularization=AbstractRegularization(),
-    )
+
+    pixelization = ag.Pixelization(mesh=AbstractMesh())
+
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[galaxy_pix], redshift=None)
     assert plane.hyper_galaxies_with_pixelization_image_list[0] is None
 
     galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=AbstractPixelization(),
-        regularization=AbstractRegularization(),
-        hyper_galaxy_image=1,
+        redshift=0.5, pixelization=pixelization, hyper_galaxy_image=1
     )
 
     plane = ag.Plane(galaxies=[galaxy_pix, ag.Galaxy(redshift=0.5)], redshift=None)
diff --git a/test_autogalaxy/plane/test_to_inversion.py b/test_autogalaxy/plane/test_to_inversion.py
index 059a8537d..b53dbf78d 100644
--- a/test_autogalaxy/plane/test_to_inversion.py
+++ b/test_autogalaxy/plane/test_to_inversion.py
@@ -4,15 +4,13 @@
 import autogalaxy as ag
 
 
-def test__lp_linear_func_galaxy_dict_from(lp_0):
+def test__lp_linear_func_list_galaxy_dict(lp_0, masked_imaging_7x7):
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5)], redshift=None)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    lp_linear_func_galaxy_dict = plane_to_inversion.lp_linear_func_galaxy_dict_from(
-        source_grid_slim=1, source_blurring_grid_slim=2, convolver=3
-    )
+    lp_linear_func_galaxy_dict = plane_to_inversion.lp_linear_func_list_galaxy_dict
 
     assert lp_linear_func_galaxy_dict == {}
 
@@ -28,11 +26,9 @@ def test__lp_linear_func_galaxy_dict_from(lp_0):
 
     plane = ag.Plane(galaxies=[galaxy_0, galaxy_1], redshift=None)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    lp_linear_func_galaxy_dict = plane_to_inversion.lp_linear_func_galaxy_dict_from(
-        source_grid_slim=1, source_blurring_grid_slim=2, convolver=3
-    )
+    lp_linear_func_galaxy_dict = plane_to_inversion.lp_linear_func_list_galaxy_dict
 
     lp_linear_func_list = list(lp_linear_func_galaxy_dict.keys())
 
@@ -40,50 +36,60 @@ def test__lp_linear_func_galaxy_dict_from(lp_0):
     assert lp_linear_func_galaxy_dict[lp_linear_func_list[1]] == galaxy_0
     assert lp_linear_func_galaxy_dict[lp_linear_func_list[2]] == galaxy_1
 
-    assert lp_linear_func_list[0].light_profile == lp_linear_0
-    assert lp_linear_func_list[1].light_profile == lp_linear_1
-    assert lp_linear_func_list[2].light_profile == lp_linear_2
+    assert lp_linear_func_list[0].light_profile_list[0] == lp_linear_0
+    assert lp_linear_func_list[1].light_profile_list[0] == lp_linear_1
+    assert lp_linear_func_list[2].light_profile_list[0] == lp_linear_2
 
+    basis = ag.lp_basis.Basis(light_profile_list=[lp_linear_0, lp_linear_1])
 
-def test__sparse_image_plane_grid_list_from(sub_grid_2d_7x7):
+    galaxy_0 = ag.Galaxy(redshift=0.5, bulge=basis)
 
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(data_pixelization_grid=[[1.0, 1.0]]),
-        regularization=ag.m.MockRegularization(),
+    plane = ag.Plane(galaxies=[galaxy_0, galaxy_1], redshift=None)
+
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
+
+    lp_linear_func_galaxy_dict = plane_to_inversion.lp_linear_func_list_galaxy_dict
+
+    lp_linear_func_list = list(lp_linear_func_galaxy_dict.keys())
+
+    assert lp_linear_func_galaxy_dict[lp_linear_func_list[0]] == galaxy_1
+    assert lp_linear_func_galaxy_dict[lp_linear_func_list[1]] == galaxy_0
+
+    assert lp_linear_func_list[0].light_profile_list[0] == lp_linear_2
+    assert lp_linear_func_list[1].light_profile_list[0] == lp_linear_0
+    assert lp_linear_func_list[1].light_profile_list[1] == lp_linear_1
+
+
+def test__sparse_image_plane_grid_list(masked_imaging_7x7):
+
+    pixelization = ag.m.MockPixelization(
+        mesh=ag.m.MockMesh(data_mesh_grid=np.array([1.0, 1.0]))
     )
 
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
+
     plane = ag.Plane(galaxies=[galaxy_pix], redshift=0.5)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    sparse_grid = plane_to_inversion.sparse_image_plane_grid_list_from(
-        grid=sub_grid_2d_7x7
-    )
+    sparse_image_plane_grid_list = plane_to_inversion.sparse_image_plane_grid_list
 
-    assert (sparse_grid == np.array([[1.0, 1.0]])).all()
+    assert (sparse_image_plane_grid_list == np.array([[1.0, 1.0]])).all()
 
     # In the ag.m.MockPixelization class the grid is returned if hyper image=None, and grid*hyper image is
     # returned otherwise.
 
     galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(
-            data_pixelization_grid=np.array([[1.0, 1.0]])
-        ),
-        regularization=ag.m.MockRegularization(),
-        hyper_galaxy_image=2,
+        redshift=0.5, pixelization=pixelization, hyper_galaxy_image=2
     )
 
     plane = ag.Plane(galaxies=[galaxy_pix], redshift=0.5)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    sparse_grid = plane_to_inversion.sparse_image_plane_grid_list_from(
-        grid=sub_grid_2d_7x7
-    )
+    sparse_image_plane_grid_list = plane_to_inversion.sparse_image_plane_grid_list
 
-    assert (sparse_grid == np.array([[2.0, 2.0]])).all()
+    assert (sparse_image_plane_grid_list == np.array([[2.0, 2.0]])).all()
 
     # No Galalxies
 
@@ -91,62 +97,51 @@ def test__sparse_image_plane_grid_list_from(sub_grid_2d_7x7):
 
     plane = ag.Plane(galaxies=[galaxy_no_pix], redshift=0.5)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    sparse_grid = plane_to_inversion.sparse_image_plane_grid_list_from(
-        grid=sub_grid_2d_7x7
-    )
+    sparse_image_plane_grid_list = plane_to_inversion.sparse_image_plane_grid_list
 
-    assert sparse_grid is None
+    assert sparse_image_plane_grid_list is None
 
 
-def test__mapper_galaxy_dict_from(sub_grid_2d_7x7):
+def test__mapper_galaxy_dict(masked_imaging_7x7):
 
-    galaxy_pix = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(
-            mapper=1, data_pixelization_grid=sub_grid_2d_7x7
-        ),
-        regularization=ag.m.MockRegularization(),
-    )
+    mesh = ag.mesh.Rectangular(shape=(3, 3))
+
+    pixelization = ag.m.MockPixelization(mesh=mesh)
+
+    galaxy_pix = ag.Galaxy(redshift=0.5, pixelization=pixelization)
     galaxy_no_pix = ag.Galaxy(redshift=0.5)
 
     plane = ag.Plane(galaxies=[galaxy_no_pix, galaxy_pix], redshift=0.5)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict_from(
-        grid=sub_grid_2d_7x7
-    )
+    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict
 
     mapper_list = list(mapper_galaxy_dict.keys())
 
-    assert mapper_list[0] == 1
+    assert mapper_list[0].pixels == 9
     assert mapper_galaxy_dict[mapper_list[0]] == galaxy_pix
 
-    galaxy_pix_2 = ag.Galaxy(
-        redshift=0.5,
-        pixelization=ag.m.MockPixelization(
-            mapper=2, data_pixelization_grid=sub_grid_2d_7x7
-        ),
-        regularization=ag.m.MockRegularization(),
-    )
+    mesh = ag.mesh.Rectangular(shape=(4, 3))
+    pixelization = ag.m.MockPixelization(mesh=mesh)
+
+    galaxy_pix_2 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
     galaxy_no_pix = ag.Galaxy(redshift=0.5)
 
     plane = ag.Plane(
         galaxies=[galaxy_no_pix, galaxy_pix, galaxy_no_pix, galaxy_pix_2], redshift=0.5
     )
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict_from(
-        grid=sub_grid_2d_7x7
-    )
+    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict
 
     mapper_list = list(mapper_galaxy_dict.keys())
 
-    assert mapper_list[0] == 1
-    assert mapper_list[1] == 2
+    assert mapper_list[0].pixels == 9
+    assert mapper_list[1].pixels == 12
 
     assert mapper_galaxy_dict[mapper_list[0]] == galaxy_pix
     assert mapper_galaxy_dict[mapper_list[1]] == galaxy_pix_2
@@ -155,52 +150,102 @@ def test__mapper_galaxy_dict_from(sub_grid_2d_7x7):
 
     plane = ag.Plane(galaxies=[galaxy_no_pix], redshift=0.5)
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
 
-    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict_from(
-        grid=sub_grid_2d_7x7
-    )
+    mapper_galaxy_dict = plane_to_inversion.mapper_galaxy_dict
 
     assert mapper_galaxy_dict == {}
 
 
+def test__regularization_list(masked_imaging_7x7):
+
+    regularization_0 = ag.reg.Constant(coefficient=1.0)
+    regularization_1 = ag.reg.ConstantSplit(coefficient=2.0)
+
+    pixelization_0 = ag.m.MockPixelization(
+        mesh=ag.mesh.Rectangular(shape=(10, 10)), regularization=regularization_0
+    )
+    pixelization_1 = ag.m.MockPixelization(
+        mesh=ag.mesh.Rectangular(shape=(8, 8)), regularization=regularization_1
+    )
+
+    galaxy_0 = ag.Galaxy(redshift=0.5, light=ag.lp_linear.EllGaussian())
+    galaxy_1 = ag.Galaxy(redshift=0.5, pixelization=pixelization_0)
+    galaxy_2 = ag.Galaxy(
+        redshift=0.5, light=ag.lp_linear.EllGaussian(), pixelization=pixelization_1
+    )
+
+    plane = ag.Plane(galaxies=[galaxy_0, galaxy_1, galaxy_2])
+
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
+
+    regularization_list = plane_to_inversion.regularization_list
+
+    assert regularization_list[0] == None
+    assert regularization_list[1] == None
+    assert regularization_list[2] == regularization_0
+    assert regularization_list[3] == regularization_1
+
+    regularization_2 = ag.reg.Constant(coefficient=3.0)
+
+    basis = ag.lp_basis.Basis(
+        light_profile_list=[ag.lp_linear.EllGaussian()], regularization=regularization_2
+    )
+
+    galaxy_3 = ag.Galaxy(redshift=0.5, bulge=basis)
+
+    plane = ag.Plane(galaxies=[galaxy_0, galaxy_1, galaxy_3])
+
+    plane_to_inversion = ag.PlaneToInversion(plane=plane, dataset=masked_imaging_7x7)
+
+    regularization_list = plane_to_inversion.regularization_list
+
+    assert regularization_list[0] == None
+    assert regularization_list[1] == regularization_2
+    assert regularization_list[2] == regularization_0
+
+
 def test__inversion_imaging_from(sub_grid_2d_7x7, masked_imaging_7x7):
 
     g_linear = ag.Galaxy(redshift=0.5, light_linear=ag.lp_linear.EllSersic())
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g_linear])
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
-
-    inversion = plane_to_inversion.inversion_imaging_from(
+    plane_to_inversion = ag.PlaneToInversion(
+        plane=plane,
         dataset=masked_imaging_7x7,
-        image=masked_imaging_7x7.image,
+        data=masked_imaging_7x7.image,
         noise_map=masked_imaging_7x7.noise_map,
         w_tilde=masked_imaging_7x7.w_tilde,
         settings_pixelization=ag.SettingsPixelization(use_border=False),
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
+    inversion = plane_to_inversion.inversion
+
     assert inversion.reconstruction[0] == pytest.approx(0.00543437, 1.0e-2)
 
-    pix = ag.pix.Rectangular(shape=(3, 3))
-    reg = ag.reg.Constant(coefficient=0.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(3, 3)),
+        regularization=ag.reg.Constant(coefficient=0.0),
+    )
 
-    g0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g0])
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
-
-    inversion = plane_to_inversion.inversion_imaging_from(
+    plane_to_inversion = ag.PlaneToInversion(
+        plane=plane,
         dataset=masked_imaging_7x7,
-        image=masked_imaging_7x7.image,
+        data=masked_imaging_7x7.image,
         noise_map=masked_imaging_7x7.noise_map,
         w_tilde=masked_imaging_7x7.w_tilde,
         settings_pixelization=ag.SettingsPixelization(use_border=False),
         settings_inversion=ag.SettingsInversion(use_w_tilde=False),
     )
 
+    inversion = plane_to_inversion.inversion
+
     assert inversion.mapped_reconstructed_image == pytest.approx(
         masked_imaging_7x7.image, 1.0e-2
     )
@@ -212,11 +257,10 @@ def test__inversion_interferometer_from(sub_grid_2d_7x7, interferometer_7):
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g_linear])
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
-
-    inversion = plane_to_inversion.inversion_interferometer_from(
+    plane_to_inversion = ag.PlaneToInversion(
+        plane=plane,
         dataset=interferometer_7,
-        visibilities=interferometer_7.visibilities,
+        data=interferometer_7.visibilities,
         noise_map=interferometer_7.noise_map,
         w_tilde=None,
         settings_pixelization=ag.SettingsPixelization(use_border=False),
@@ -225,22 +269,25 @@ def test__inversion_interferometer_from(sub_grid_2d_7x7, interferometer_7):
         ),
     )
 
+    inversion = plane_to_inversion.inversion
+
     assert inversion.reconstruction[0] == pytest.approx(0.0012073, 1.0e-2)
 
     interferometer_7.data = ag.Visibilities.ones(shape_slim=(7,))
 
-    pix = ag.pix.Rectangular(shape=(7, 7))
-    reg = ag.reg.Constant(coefficient=0.0)
+    pixelization = ag.Pixelization(
+        mesh=ag.mesh.Rectangular(shape=(7, 7)),
+        regularization=ag.reg.Constant(coefficient=0.0),
+    )
 
-    g0 = ag.Galaxy(redshift=0.5, pixelization=pix, regularization=reg)
+    g0 = ag.Galaxy(redshift=0.5, pixelization=pixelization)
 
     plane = ag.Plane(galaxies=[ag.Galaxy(redshift=0.5), g0])
 
-    plane_to_inversion = ag.PlaneToInversion(plane=plane)
-
-    inversion = plane_to_inversion.inversion_interferometer_from(
+    plane_to_inversion = ag.PlaneToInversion(
+        plane=plane,
         dataset=interferometer_7,
-        visibilities=interferometer_7.visibilities,
+        data=interferometer_7.visibilities,
         noise_map=interferometer_7.noise_map,
         w_tilde=None,
         settings_pixelization=ag.SettingsPixelization(use_border=False),
@@ -249,6 +296,8 @@ def test__inversion_interferometer_from(sub_grid_2d_7x7, interferometer_7):
         ),
     )
 
+    inversion = plane_to_inversion.inversion
+
     assert inversion.mapped_reconstructed_data.real == pytest.approx(
         interferometer_7.visibilities.real, 1.0e-2
     )
diff --git a/test_autogalaxy/plot/mat_wrap/config/visualize/include.ini b/test_autogalaxy/plot/mat_wrap/config/visualize/include.ini
index 8318adccc..d09237994 100644
--- a/test_autogalaxy/plot/mat_wrap/config/visualize/include.ini
+++ b/test_autogalaxy/plot/mat_wrap/config/visualize/include.ini
@@ -8,8 +8,8 @@ caustics=False
 multiple_images=False
 light_profile_centres=True
 mass_profile_centres=True
-mapper_data_pixelization_grid=True
-mapper_source_pixelization_grid=True
+mapper_data_mesh_grid=True
+mapper_source_mesh_grid=False
 mapper_source_grid_slim=False
 parallel_overscan=True
 serial_prescan=True
diff --git a/test_autogalaxy/profiles/config/summary/priors/test_analysis.json b/test_autogalaxy/profiles/config/summary/priors/test_analysis.json
index 16bac7a1e..32fac4991 100644
--- a/test_autogalaxy/profiles/config/summary/priors/test_analysis.json
+++ b/test_autogalaxy/profiles/config/summary/priors/test_analysis.json
@@ -1,5 +1,5 @@
 {
-    "aa.m.MockPixelization": {
+    "aa.m.MockMesh": {
         "number_clusters": {
             "type": "Uniform",
             "lower_limit": 1.0,
diff --git a/test_autogalaxy/profiles/light_profiles/test_basis.py b/test_autogalaxy/profiles/light_profiles/test_basis.py
new file mode 100644
index 000000000..67281882a
--- /dev/null
+++ b/test_autogalaxy/profiles/light_profiles/test_basis.py
@@ -0,0 +1,57 @@
+import numpy as np
+import pytest
+
+import autogalaxy as ag
+
+
+def test__image_2d_from__does_not_include_linear_light_profiles(sub_grid_2d_7x7):
+
+    lp = ag.lp.EllSersic(intensity=0.1)
+
+    lp_linear = ag.lp_linear.EllSersic(effective_radius=2.0, sersic_index=2.0)
+
+    lp_image = lp.image_2d_from(grid=sub_grid_2d_7x7)
+
+    basis = ag.lp_basis.Basis(light_profile_list=[lp, lp_linear])
+
+    image = basis.image_2d_from(grid=sub_grid_2d_7x7)
+
+    assert (image == lp_image).all()
+
+
+def test__image_2d_from__operated_only_input(sub_grid_2d_7x7, lp_0, lp_operated_0):
+
+    image_2d_not_operated = lp_0.image_2d_from(grid=sub_grid_2d_7x7)
+    image_2d_operated = lp_operated_0.image_2d_from(grid=sub_grid_2d_7x7)
+
+    basis = ag.lp_basis.Basis(light_profile_list=[lp_0, lp_operated_0])
+
+    image_2d = basis.image_2d_from(grid=sub_grid_2d_7x7, operated_only=False)
+    assert (image_2d == image_2d_not_operated).all()
+
+    image_2d = basis.image_2d_from(grid=sub_grid_2d_7x7, operated_only=True)
+    assert (image_2d == image_2d_operated).all()
+
+    image_2d = basis.image_2d_from(grid=sub_grid_2d_7x7, operated_only=None)
+    assert (image_2d == image_2d_not_operated + image_2d_operated).all()
+
+
+def test__image_2d_list_from__operated_only_input(sub_grid_2d_7x7, lp_0, lp_operated_0):
+
+    image_2d_not_operated = lp_0.image_2d_from(grid=sub_grid_2d_7x7)
+    image_2d_operated = lp_operated_0.image_2d_from(grid=sub_grid_2d_7x7)
+
+    basis = ag.lp_basis.Basis(light_profile_list=[lp_0, lp_operated_0])
+
+    image_2d_list = basis.image_2d_list_from(grid=sub_grid_2d_7x7, operated_only=False)
+    assert (image_2d_list[0] == image_2d_not_operated).all()
+    assert (image_2d_list[1] == np.zeros((36))).all()
+
+    image_2d_list = basis.image_2d_list_from(grid=sub_grid_2d_7x7, operated_only=True)
+    assert (image_2d_list[0] == np.zeros((36))).all()
+    assert (image_2d_list[1] == image_2d_operated).all()
+
+    image_2d_list = basis.image_2d_list_from(grid=sub_grid_2d_7x7, operated_only=None)
+    assert (
+        image_2d_list[0] + image_2d_list[1] == image_2d_not_operated + image_2d_operated
+    ).all()
diff --git a/test_autogalaxy/profiles/light_profiles/test_light_profiles_linear.py b/test_autogalaxy/profiles/light_profiles/test_light_profiles_linear.py
index 9e2ff9d9a..e789b85bb 100644
--- a/test_autogalaxy/profiles/light_profiles/test_light_profiles_linear.py
+++ b/test_autogalaxy/profiles/light_profiles/test_light_profiles_linear.py
@@ -1,6 +1,56 @@
+import numpy as np
+import pytest
+
 import autogalaxy as ag
 
 from autogalaxy.profiles.light_profiles.light_profiles_linear import LightProfileLinear
+from autogalaxy.profiles.light_profiles.light_profiles_linear import (
+    LightProfileLinearObjFuncList,
+)
+
+
+def test__mapping_matrix_from(sub_grid_2d_7x7, blurring_grid_2d_7x7, convolver_7x7):
+
+    lp_0 = ag.lp_linear.EllSersic(effective_radius=1.0)
+    lp_1 = ag.lp_linear.EllSersic(effective_radius=2.0)
+
+    lp_linear_obj_func_list = LightProfileLinearObjFuncList(
+        grid=sub_grid_2d_7x7,
+        blurring_grid=blurring_grid_2d_7x7,
+        convolver=convolver_7x7,
+        light_profile_list=[lp_0, lp_1],
+    )
+
+    assert lp_linear_obj_func_list.pixels == 2
+
+    lp_0_image = lp_0.image_2d_from(grid=sub_grid_2d_7x7)
+    lp_1_image = lp_1.image_2d_from(grid=sub_grid_2d_7x7)
+
+    assert lp_linear_obj_func_list.mapping_matrix[:, 0] == pytest.approx(
+        lp_0_image.binned, 1.0e-4
+    )
+    assert lp_linear_obj_func_list.mapping_matrix[:, 1] == pytest.approx(
+        lp_1_image.binned, 1.0e-4
+    )
+
+    lp_0_blurred_image = lp_0.blurred_image_2d_from(
+        grid=sub_grid_2d_7x7,
+        blurring_grid=blurring_grid_2d_7x7,
+        convolver=convolver_7x7,
+    )
+
+    lp_1_blurred_image = lp_1.blurred_image_2d_from(
+        grid=sub_grid_2d_7x7,
+        blurring_grid=blurring_grid_2d_7x7,
+        convolver=convolver_7x7,
+    )
+
+    assert lp_linear_obj_func_list.operated_mapping_matrix_override[
+        :, 0
+    ] == pytest.approx(lp_0_blurred_image, 1.0e-4)
+    assert lp_linear_obj_func_list.operated_mapping_matrix_override[
+        :, 1
+    ] == pytest.approx(lp_1_blurred_image, 1.0e-4)
 
 
 def test__lp_from():
diff --git a/test_autogalaxy/profiles/mass_profiles/config/summary/priors/test_analysis.json b/test_autogalaxy/profiles/mass_profiles/config/summary/priors/test_analysis.json
index 16bac7a1e..32fac4991 100644
--- a/test_autogalaxy/profiles/mass_profiles/config/summary/priors/test_analysis.json
+++ b/test_autogalaxy/profiles/mass_profiles/config/summary/priors/test_analysis.json
@@ -1,5 +1,5 @@
 {
-    "aa.m.MockPixelization": {
+    "aa.m.MockMesh": {
         "number_clusters": {
             "type": "Uniform",
             "lower_limit": 1.0,