diff --git a/projects/README.rst b/projects/README.rst
new file mode 100644
index 0000000..3226487
--- /dev/null
+++ b/projects/README.rst
@@ -0,0 +1,7 @@
+The ``projects`` contains example projects which use **PyAutoFit** to perform modeling.
+
+They show how one would build a structure a project around **PyAutoFit** to best use its tools for model-composition, model-fitting, visualization, etc.
+
+They also demonstrate advanced functionality using real-world examples.
+
+- ``cosmology``: Example **PyAutoFit** project using a model-fitting example from Astronomy. This highlights **multi-level models**.
\ No newline at end of file
diff --git a/projects/__init__.py b/projects/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/projects/cosmology/__init__.py b/projects/cosmology/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/projects/cosmology/config/.gitignore b/projects/cosmology/config/.gitignore
new file mode 100644
index 0000000..b722e9e
--- /dev/null
+++ b/projects/cosmology/config/.gitignore
@@ -0,0 +1 @@
+!.gitignore
\ No newline at end of file
diff --git a/projects/cosmology/config/README.rst b/projects/cosmology/config/README.rst
new file mode 100644
index 0000000..7c08fcf
--- /dev/null
+++ b/projects/cosmology/config/README.rst
@@ -0,0 +1,15 @@
+The ``config`` folder contains configuration files which customize default **PyAutoLens**.
+
+Folders
+-------
+
+- ``priors``: Configs defining default priors assumed on every model component and set of parameters.
+
+Files
+-----
+
+- ``general.yaml``: Customizes general **PyAutoLens** settings.
+- ``non-linear.yaml``: Configs for default non-linear search (e.g. MCMC, nested sampling) settings.
+- ``logging.yaml``: Customizes the logging behaviour of **PyAutoLens**.
+- ``visualize.yaml``: Configs defining what images are output by a lens model fit.
+- ``notation.yaml``: Configs defining labels and formatting of model parameters when used for visualization.
diff --git a/projects/cosmology/config/general.yaml b/projects/cosmology/config/general.yaml
new file mode 100644
index 0000000..ca17db7
--- /dev/null
+++ b/projects/cosmology/config/general.yaml
@@ -0,0 +1,33 @@
+updates:
+  iterations_per_quick_update: 1e99 # Non-linear search iterations between every quick update, which just displays the maximum likelihood model fit.
+  iterations_per_full_update: 1e99  # Non-linear search iterations between every full update, which outputs all visuals and result fits (e.g. model.result, search.summary), this exits the search and can be slow.  
+hpc:
+  hpc_mode: false                   # If True, use HPC mode, which disables GUI visualization, logging to screen and other settings which are not suited to running on a super computer.
+  iterations_per_quick_update: 10000 #  Non-linear search iterations between every quick update, which just displays the maximum likelihood model fit.  
+  iterations_per_full_update: 1e99  # Non-linear search iterations between every full update, which outputs all visuals and result fits (e.g. model.result, search.summary), this exits the search and can be slow.
+inversion:
+  check_reconstruction: true        # If True, the inversion's reconstruction is checked to ensure the solution of a meshs's mapper is not an invalid solution where the values are all the same.
+  reconstruction_vmax_factor: 0.5   # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.  
+output:
+  force_pickle_overwrite: false     #   force_pickle_overwrite: false     # If True, pickle files output by a search (e.g. samples.pickle) are recreated when a new model-fit is performed.
+  force_visualize_overwrite: false # If True, visualization images output by a search (e.g. subplots of the fit) are recreated when a new model-fit is performed.
+  info_whitespace_length: 80        # Length of whitespace between the parameter names and values in the model.info / result.info
+  log_level: INFO                   # The level of information output by logging.
+  log_to_file: false                # If True, outputs the non-linear search log to a file (and not printed to screen).
+  log_file: output.log              # The name of the file the logged output is written to (in the non-linear search output folder)
+  model_results_decimal_places: 3   # Number of decimal places estimated parameter values / errors are output in model.results.
+  remove_files: false               # If True, all output files of a non-linear search (e.g. samples, visualization, etc.) are deleted once the model-fit has completed, such that only the .zip file remains.
+  search_internal : false    # If True, the search internal folder which contains a saved state of the non-linear search is delete, saving on hard-disk space.
+  samples_to_csv: true              # If True, non-linear search samples are written to a .csv file.
+  unconverged_sample_size : 100     # If outputting results of an unconverged search, the number of samples used to estimate the median PDF values and errors.
+parallel:
+  warn_environment_variables: true  # If True, a warning is displayed when the search's number of CPU > 1 and enviromment variables related to threading are also > 1.
+
+profiling:
+  parallel_profile: false           # If True, the parallelization of the fit is profiled outputting a cPython graph.
+  should_profile: false             # If True, the ``profile_log_likelihood_function()`` function of an analysis class is called throughout a model-fit, profiling run times.
+  repeats: 1                        # The number of repeat function calls used to measure run-times when profiling.
+test:
+  check_preloads: false
+  exception_override: false
+  parallel_profile: false
diff --git a/projects/cosmology/config/logging.yaml b/projects/cosmology/config/logging.yaml
new file mode 100644
index 0000000..273f702
--- /dev/null
+++ b/projects/cosmology/config/logging.yaml
@@ -0,0 +1,20 @@
+version: 1
+disable_existing_loggers: false
+handlers:
+  console:
+    class: logging.StreamHandler
+    level: INFO
+    stream: ext://sys.stdout
+    formatter: formatter
+  file:
+    class: logging.FileHandler
+    level: INFO
+    filename: root.log
+    formatter: formatter
+root:
+  level: INFO
+  handlers:
+  - console
+formatters:
+  formatter:
+    format: '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
diff --git a/projects/cosmology/config/non_linear/GridSearch.yaml b/projects/cosmology/config/non_linear/GridSearch.yaml
new file mode 100644
index 0000000..af6daba
--- /dev/null
+++ b/projects/cosmology/config/non_linear/GridSearch.yaml
@@ -0,0 +1,5 @@
+# The settings of a parallelized grid search of non-linear searches.
+
+parallel:
+  number_of_cores: 3 # The number of cores the search is parallelized over by default, using Python multiprocessing.
+  step_size: 0.1     # The default step size of each grid search parameter, in terms of unit values of the priors.
\ No newline at end of file
diff --git a/projects/cosmology/config/non_linear/README.rst b/projects/cosmology/config/non_linear/README.rst
new file mode 100644
index 0000000..9432dd8
--- /dev/null
+++ b/projects/cosmology/config/non_linear/README.rst
@@ -0,0 +1,9 @@
+The ``non_linear`` folder contains configuration files which customize the default behaviour of non-linear searches in
+**PyAutoLens**.
+
+Files
+-----
+
+- ``mcmc.yaml``: Settings default behaviour of MCMC non-linear searches (e.g. Emcee).
+- ``nest.yaml``: Settings default behaviour of nested sampler non-linear searches (e.g. Dynesty).
+- ``mle.yaml``: Settings default behaviour of maximum likelihood estimator (mle) searches (e.g. PySwarms).
\ No newline at end of file
diff --git a/projects/cosmology/config/notation.yaml b/projects/cosmology/config/notation.yaml
new file mode 100644
index 0000000..8c525ca
--- /dev/null
+++ b/projects/cosmology/config/notation.yaml
@@ -0,0 +1,41 @@
+# The notation configs define the labels of every model parameter and its derived quantities, which are used when
+# visualizing results (for example labeling the axis of the PDF triangle plots output by a non-linear search).
+
+
+# label: The label given to the each parameter, for plots like PDF corner plots.
+
+# For example, if `centre=x`, the plot axis will be labeled 'x'.
+
+
+# superscript: the superscript used on certain plots that show the results of different model-components.
+
+# For example, if `Gaussian=g`, plots where the parameters of the Gaussian model-component have superscript `g`.
+
+label:
+  label:
+    sigma: \sigma
+    centre: x
+    intensity: norm
+    parameter0: a
+    parameter1: b
+    parameter2: c
+    rate: \lambda
+  superscript:
+    Exponential: e
+    Gaussian: g
+    ModelComponent0: M0
+    ModelComponent1: M1
+
+# label_format: The format certain parameters are output as in output files like the `model.results` file.
+
+# For example, if `centre={:.2d}`, the format of the centre parameter in results files will use this Python format.
+
+label_format:
+  format:
+    sigma: '{:.2f}'
+    centre: '{:.2f}'
+    intensity: '{:.2f}'
+    parameter0: '{:.2f}'
+    parameter1: '{:.2f}'
+    parameter2: '{:.2f}'
+    rate: '{:.2f}'
diff --git a/projects/cosmology/config/priors/galaxy.yaml b/projects/cosmology/config/priors/galaxy.yaml
new file mode 100644
index 0000000..1938e8a
--- /dev/null
+++ b/projects/cosmology/config/priors/galaxy.yaml
@@ -0,0 +1,11 @@
+Redshift:
+  redshift:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 3.0
+    width_modifier:
+      type: Absolute
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
\ No newline at end of file
diff --git a/projects/cosmology/config/priors/light_profiles.yaml b/projects/cosmology/config/priors/light_profiles.yaml
new file mode 100644
index 0000000..4e10fb9
--- /dev/null
+++ b/projects/cosmology/config/priors/light_profiles.yaml
@@ -0,0 +1,61 @@
+LightDevVaucouleurs:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.3
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.3
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  axis_ratio:
+    type: Uniform
+    lower_limit: 0.2
+    upper_limit: 1.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
+  angle:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 180.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
+  effective_radius:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 30.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
+  intensity:
+    type: LogUniform
+    lower_limit: 1.0e-06
+    upper_limit: 1000000.0
+    width_modifier:
+      type: Relative
+      value: 0.5
+    limits:
+      lower: 0.0
+      upper: inf
diff --git a/projects/cosmology/config/priors/mass_profiles.yaml b/projects/cosmology/config/priors/mass_profiles.yaml
new file mode 100644
index 0000000..9c01a33
--- /dev/null
+++ b/projects/cosmology/config/priors/mass_profiles.yaml
@@ -0,0 +1,51 @@
+MassIsothermal:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.3
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.3
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  axis_ratio:
+    type: Uniform
+    lower_limit: 0.2
+    upper_limit: 1.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
+  angle:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 180.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
+  mass:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 4.0
+    width_modifier:
+      type: Relative
+      value: 1.0
+    limits:
+      lower: 0.0
+      upper: inf
diff --git a/projects/cosmology/config/visualize.yaml b/projects/cosmology/config/visualize.yaml
new file mode 100644
index 0000000..64cfa0d
--- /dev/null
+++ b/projects/cosmology/config/visualize.yaml
@@ -0,0 +1,28 @@
+general:
+  general:
+    backend: default         # The matploblib backend used for visualization. `default` uses the system default, can specifiy specific backend (e.g. TKAgg, Qt5Agg, WXAgg).
+plots_search:
+  dynesty:
+    corner: true             # Output corner figure (using anestetic) during a non-linear search fit?
+    cornerpoints: false      # Output Dynesty cornerpoints figure during a non-linear search fit?
+    runplot: true            # Output Dynesty runplot figure during a non-linear search fit?
+    traceplot: true          # Output Dynesty traceplot figure during a non-linear search fit?
+  emcee:
+    corner: true             # Output corner figure (using corner.py) during a non-linear search fit?
+    likelihood_series: true  # Output likelihood series figure during a non-linear search fit?
+    time_series: true        # Output time series figure during a non-linear search fit?
+    trajectories: true       # Output trajectories figure during a non-linear search fit?
+  pyswarms:
+    contour: true            # Output contour figure during a non-linear search fit?
+    cost_history: true       # Output cost_history figure during a non-linear search fit?
+    time_series: true        # Output time_series figure during a non-linear search fit?
+    trajectories: true       # Output trajectories figure during a non-linear search fit?
+  ultranest:
+    corner: true             # Output corner figure (using anestetic) during a non-linear search fit?
+    runplot: true            # Output Ultranest runplot figure during a non-linear search fit?
+    traceplot: true          # Output Ultranest traceplot figure during a non-linear search fit?
+  zeus:
+    corner: true             # Output corner figure (using corner.py) during a non-linear search fit?
+    likelihood_series: true  # Output likelihood series figure during a non-linear search fit?
+    time_series: true        # Output time series figure during a non-linear search fit?
+    trajectories: true       # Output trajectories figure during a non-linear search fit?
\ No newline at end of file
diff --git a/projects/cosmology/dataset/data.npy b/projects/cosmology/dataset/data.npy
new file mode 100644
index 0000000..229ad12
Binary files /dev/null and b/projects/cosmology/dataset/data.npy differ
diff --git a/projects/cosmology/dataset/grid.npy b/projects/cosmology/dataset/grid.npy
new file mode 100644
index 0000000..8b62f92
Binary files /dev/null and b/projects/cosmology/dataset/grid.npy differ
diff --git a/projects/cosmology/dataset/noise_map.npy b/projects/cosmology/dataset/noise_map.npy
new file mode 100644
index 0000000..a136534
Binary files /dev/null and b/projects/cosmology/dataset/noise_map.npy differ
diff --git a/projects/cosmology/dataset/psf.npy b/projects/cosmology/dataset/psf.npy
new file mode 100644
index 0000000..3274867
Binary files /dev/null and b/projects/cosmology/dataset/psf.npy differ
diff --git a/projects/cosmology/example_1_intro.ipynb b/projects/cosmology/example_1_intro.ipynb
new file mode 100644
index 0000000..eec0de0
--- /dev/null
+++ b/projects/cosmology/example_1_intro.ipynb
@@ -0,0 +1,583 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Project: Cosmology\n",
+        "==================\n",
+        "\n",
+        "This project uses the astrophysical phenomena of Strong Gravitational Lensing to illustrate multi-level model\n",
+        "composition and fitting with **PyAutoFit**.\n",
+        "\n",
+        "A strong gravitational lens is a system where two (or more) galaxies align perfectly down our line of sight from Earth\n",
+        "such that the foreground galaxy's mass deflects the light of a background source galaxy(s).\n",
+        "\n",
+        "When the alignment is just right and the lens is massive enough, the background source galaxy appears multiple\n",
+        "times. The schematic below shows such a system, where light-rays from the source are deflected around the lens galaxy\n",
+        "to the observer following multiple distinct paths.\n",
+        "\n",
+        "![Schematic of Gravitational Lensing](https://raw.githubusercontent.com/Jammy2211/PyAutoLens/main/docs/overview/images/overview_1_lensing/schematic.jpg)\n",
+        "**Credit: F. Courbin, S. G. Djorgovski, G. Meylan, et al., Caltech / EPFL / WMKO**\n",
+        "https://www.cosmology.caltech.edu/~george/qsolens/\n",
+        "\n",
+        "As an observer, we don't see the source's true appearance (e.g. the red round blob of light). We only observe its\n",
+        "light after it has been deflected and lensed by the foreground galaxies (e.g. as the two distinct red multiple images\n",
+        " in the image on the left). We also observe the emission of the foreground galaxy (in blue).\n",
+        "\n",
+        "You can read more about gravitational lensing as the following link:\n",
+        "\n",
+        "https://en.wikipedia.org/wiki/Gravitational_lens\n",
+        "\n",
+        "__PyAutoLens__\n",
+        "\n",
+        "Strong gravitational lensing is the original science case that sparked the development of **PyAutoFit**, which is\n",
+        "a spin off of our astronomy software **PyAutoLens** `https://github.com/Jammy2211/PyAutoLens`.\n",
+        "\n",
+        "We'll use **PyAutoLens** to illustrate how the tools we developed with **PyAutoFit** allowed us to\n",
+        "ensure **PyAutoLens**'s model fitting tools were extensible, easy to maintain and enabled intuitive model composition.\n",
+        "\n",
+        "__Multi-Level Models__\n",
+        "\n",
+        "Strong lensing is a great case study for using **PyAutoFit**, due to the multi-component nature of how one composes\n",
+        "a strong lens model. A strong lens model consists of light and mass models of each galaxy in the lens system, where\n",
+        "each galaxy is a model in itself. The galaxies are combined into one overall \"lens model\", which in later tutorials\n",
+        "we will show may also have a Cosmological model.\n",
+        "\n",
+        "This example project uses **PyAutoFit** to compose and fit models of a strong lens, in particular highlighting\n",
+        "**PyAutoFits** multi-level model composition.\n",
+        "\n",
+        "__Strong Lens Modeling__\n",
+        "\n",
+        "The models are fitted to Hubble Space Telescope imaging of a real strong lens system and will allow us to come up\n",
+        "with a description of how light is deflected on its path through the Universe.\n",
+        "\n",
+        "This project consists of two example scripts / notebooks:\n",
+        "\n",
+        " 1) `example_1_intro`: An introduction to strong lensing, and the various parts of the project's source code that are\n",
+        "    used to represent a strong lens galaxy.\n",
+        "\n",
+        " 2) `example_2_multi_level_model`: Using **PyAutoFit** to model a strong lens, with a strong emphasis on the\n",
+        "    multi-level model API.\n",
+        "\n",
+        "__This Example__\n",
+        "\n",
+        "This introduction primarily focuses on what strong lensing is, how we define the individual model-components and fit\n",
+        "a strong lens model to data. It does not make much use of **PyAutoFit**, but it does provide a clear understanding of\n",
+        "the model so that **PyAutoFit**'s use in example 2 is clear.\n",
+        "\n",
+        "Note that import `import src as cosmo`. The package `src` contains all the code we need for this example Cosmology\n",
+        "use case, and can be thought of as the source-code you would write to perform model-fitting via **PyAutoFit** for your\n",
+        "problem of interest.\n",
+        "\n",
+        "The module `src/__init__.py` performs a series of imports that are used throughout this lecture to provide convenient\n",
+        "access to different parts of the source code."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "# %matplotlib inline\n",
+        "# from pyprojroot import here\n",
+        "# workspace_path = str(here())\n",
+        "# %cd $workspace_path\n",
+        "# print(f\"Working Directory has been set to `{workspace_path}`\")\n",
+        "\n",
+        "import src as cosmo\n",
+        "import matplotlib.pyplot as plt\n",
+        "import numpy as np\n",
+        "from scipy import signal\n",
+        "from os import path"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Plot__\n",
+        "\n",
+        "We will plot a lot of arrays of 2D data and grids of 2D coordinates in this example, so lets make a convenience \n",
+        "functions."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "\n",
+        "def plot_array(array, title=None, norm=None):\n",
+        "    plt.imshow(array, norm=norm)\n",
+        "    plt.colorbar()\n",
+        "    plt.title(title)\n",
+        "    plt.show()\n",
+        "    plt.close()\n",
+        "\n",
+        "\n",
+        "def plot_grid(grid, title=None):\n",
+        "    plt.scatter(x=grid[:, :, 0], y=grid[:, :, 1], s=1)\n",
+        "    plt.title(title)\n",
+        "    plt.show()\n",
+        "    plt.close()\n"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Data__\n",
+        "\n",
+        "First, lets load and plot Hubble Space Telescope imaging data of the strong gravitational lens called SDSSJ2303+1422, \n",
+        "where this data includes:\n",
+        " \n",
+        " 1) The image of the strong lens, which is the data we'll fit.\n",
+        " 2) The noise in every pixel of this image, which will be used when evaluating the log likelihood.\n",
+        "\n",
+        "__Masking__\n",
+        "\n",
+        "When fitting 2D imaging data, it is common to apply a mask which removes regions of the image that are not relevant to\n",
+        "the model fitting.\n",
+        "\n",
+        "For example, when fitting the strong lens, we remove the edges of the image where the lens and source galaxy's light is \n",
+        "not visible.\n",
+        "\n",
+        "In the strong lens image and noise map below, you can see this has already been performed, with the edge regions\n",
+        "blank."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "dataset_path = path.join(\"projects\", \"cosmology\", \"dataset\")\n",
+        "\n",
+        "data = np.load(file=path.join(dataset_path, \"data.npy\"))\n",
+        "plot_array(array=data, title=\"Image of Strong Lens SDSSJ2303+1422\")\n",
+        "\n",
+        "noise_map = np.load(file=path.join(dataset_path, \"noise_map.npy\"))\n",
+        "plot_array(array=noise_map, title=\"Noise Map of Strong Lens SDSSJ2303+1422\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "In the image of the strong lens two distinct objects can be seen:\n",
+        "\n",
+        " 1) A central blob of light, corresponding to the foreground lens galaxy whose mass is responsible for deflecting light.\n",
+        " 2) Two faint arcs of light in the bakcground, which is the lensed background source.\n",
+        " \n",
+        "__PSF__\n",
+        "\n",
+        "Another component of imaging data is the Point Spread Function (PSF), which describes how the light of the galaxies\n",
+        "are blurred when they enter the Huble Space Telescope's. \n",
+        "\n",
+        "This is because diffraction occurs when the light enters HST's optics, causing the light to smear out. The PSF is\n",
+        "a two dimensional array that describes this blurring via a 2D convolution kernel.\n",
+        "\n",
+        "When fitting the data below and in the `log_likelihood_function`, you'll see that the PSF is used when creating the \n",
+        "model data. This is an example of how an `Analysis` class may be extended to include additional steps in the model\n",
+        "fitting procedure."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "psf = np.load(file=path.join(dataset_path, \"psf.npy\"))\n",
+        "plot_array(array=psf, title=\"Point Spread Function of Strong Lens SDSSJ2303+1422\")\n"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Grid__\n",
+        " \n",
+        "To perform strong lensing, we need a grid of (x,y) coordinates which we map throughout the Universe as if their path\n",
+        "is deflected. \n",
+        "\n",
+        "For this, we create a simple 2D grid of coordinates below where the origin is (0.0, 0.0) and the size of\n",
+        "a pixel is 0.05, which corresponds to the resolution of our image `data`. \n",
+        "\n",
+        "This grid only contains (y,x) coordinates within the cricular mask that was applied to the data, as we only need to\n",
+        "perform ray-tracing within this region."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "grid = np.load(file=path.join(dataset_path, \"grid.npy\"))\n",
+        "\n",
+        "plot_grid(\n",
+        "    grid=grid,\n",
+        "    title=\"Cartesian grid of (x,y) coordinates aligned with strong lens dataset\",\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Light Profiles__\n",
+        "\n",
+        "Our model of a strong lens must include a description of the light of each galaxy, which we call a \"light profile\".\n",
+        "In the source-code of this example project, specifically the module `src/light_profiles.py` you will see there\n",
+        "are two light profile classes named `LightDeVaucouleurs` and `LightExponential`.\n",
+        "\n",
+        "These Python classes are the model components we will use to represent each galaxy's light and they behave analogous \n",
+        "to the `Gaussian` class seen in other tutorials. The input parameters of their `__init__` constructor (e.g. `centre`, \n",
+        "`axis_ratio`, `angle`) are their model parameters that may be fitted for by a non-linear search.\n",
+        "\n",
+        "These classes also contain functions which create an image from the light profile if an input grid of (x,y) 2D \n",
+        "coordinates are input, which we use below to create an image of a light profile."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "light_profile = cosmo.lp.LightExponential(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, intensity=1.0, effective_radius=2.0\n",
+        ")\n",
+        "light_image = light_profile.image_from_grid(grid=grid)\n",
+        "\n",
+        "plot_array(array=light_image, title=\"Image of an Exponential light profile.\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Mass Profiles__\n",
+        "\n",
+        "Our model also includes the mass of the foreground lens galaxy, called a 'mass profile'. In the source-code of the \n",
+        "example project, specifically the module `src/mass_profiles.py` you will see there is a mass profile class named \n",
+        "`MassIsothermal`. Like the light profile, this will be a model-component **PyAutoFit** fits via a non-linear search.\n",
+        "\n",
+        "The class also contains functions which create the \"deflections angles\", which describe the angles by which light is \n",
+        "deflected when it passes the mass of the foreground lens galaxy. These are subtracted from the (y,x) grid above to\n",
+        "determine the original coordinates of the source galaxy before lensing.\n",
+        "\n",
+        "A higher mass galaxy, which bends light more, will have higher values of the deflection angles plotted below:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "mass_profile = cosmo.mp.MassIsothermal(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.5\n",
+        ")\n",
+        "mass_deflections = mass_profile.deflections_from_grid(grid=grid)\n",
+        "\n",
+        "plot_array(\n",
+        "    array=mass_deflections[:, :, 0],\n",
+        "    title=\"X-component of the deflection angles of a Isothermal mass profile.\",\n",
+        ")\n",
+        "plot_array(\n",
+        "    array=mass_deflections[:, :, 1],\n",
+        "    title=\"Y-component of the deflection angles of a Isothermal mass profile.\",\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Ray Tracing__\n",
+        "\n",
+        "The deflection angles describe how our (x,y) grid of coordinates are deflected by the mass of the foreground galaxy.\n",
+        "\n",
+        "We can therefore ray-trace the grid aligned with SDSSJ2303+1422 using the mass profile above and plot a grid of\n",
+        "coordinates in the reference frame of before their light is gravitationally lensed:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "traced_grid = grid - mass_deflections\n",
+        "\n",
+        "plot_grid(grid=traced_grid, title=\"Cartesian grid of (x,y) traced coordinates.\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "By inputting this traced grid of (x,y) coordinates into our light profile, we can create an image of the galaxy as if\n",
+        "it were gravitationally lensed by the mass profile."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "traced_light_image = light_profile.image_from_grid(grid=traced_grid)\n",
+        "\n",
+        "plot_array(\n",
+        "    array=traced_light_image,\n",
+        "    title=\"Image of a gravitationally lensed Exponential light profile.\",\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Galaxy__\n",
+        "\n",
+        "In the `src/galaxy.py` module we define the `Galaxy` class, which is a collection of light and mass profiles at an \n",
+        "input redshift. For strong lens modeling, we have to use `Galaxy` objects, as the redshifts define how ray-tracing is\n",
+        "performed.\n",
+        "\n",
+        "We create two instances of the `Galaxy` class, representing the lens and source galaxies in a strong lens system."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "light_profile = cosmo.lp.LightDeVaucouleurs(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.9, angle=45.0, intensity=0.1, effective_radius=1.0\n",
+        ")\n",
+        "mass_profile = cosmo.mp.MassIsothermal(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.8\n",
+        ")\n",
+        "lens_galaxy = cosmo.Galaxy(\n",
+        "    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]\n",
+        ")\n",
+        "\n",
+        "light_profile = cosmo.lp.LightExponential(\n",
+        "    centre=(0.1, 0.1), axis_ratio=0.5, angle=80.0, intensity=1.0, effective_radius=5.0\n",
+        ")\n",
+        "source_galaxy = cosmo.Galaxy(\n",
+        "    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "A galaxy's image is the sum of its light profile images, and its deflection angles are the sum of its mass profile\n",
+        "deflection angles.\n",
+        "\n",
+        "To illustrate this, lets plot the lens galaxy's light profile image:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "galaxy_image = lens_galaxy.image_from_grid(grid=grid)\n",
+        "\n",
+        "plot_array(array=galaxy_image, title=\"Image of the Lens Galaxy.\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Data Fitting__\n",
+        "\n",
+        "We can create an overall image of the strong lens by:\n",
+        "\n",
+        " 1) Creating an image of the lens galaxy.\n",
+        " 2) Computing the deflection angles of the lens galaxy.\n",
+        " 3) Ray-tracing light to the source galaxy reference frame and using its light profile to make its image."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "lens_image = lens_galaxy.image_from_grid(grid=grid)\n",
+        "lens_deflections = lens_galaxy.deflections_from_grid(grid=grid)\n",
+        "\n",
+        "traced_grid = grid - lens_deflections\n",
+        "\n",
+        "source_image = source_galaxy.image_from_grid(grid=traced_grid)\n",
+        "\n",
+        "# The grid has zeros at its edges, which produce nans in the model image.\n",
+        "# These lead to an ill-defined log likelihood, so we set them to zero.\n",
+        "overall_image = np.nan_to_num(overall_image)\n",
+        "\n",
+        "plot_array(array=overall_image, title=\"Image of the overall Strong Lens System.\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Model Data__\n",
+        "\n",
+        "To produce the `model_data`, we now convolution the overall image with the Point Spread Function (PSF) of our\n",
+        "observations. This blurs the image to simulate the telescope optics and pixelization used to observe the image."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "model_data = signal.convolve2d(overall_image, psf, mode=\"same\")\n",
+        "\n",
+        "\n",
+        "plot_array(array=model_data, title=\"Image of the overall Strong Lens System.\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "By subtracting this model image from the data, we can create a 2D residual map. This is equivalent to the residual maps\n",
+        "we made in the 1D Gaussian examples, except for 2D imaging data.\n",
+        "\n",
+        "Clearly, the random lens model we used in this example does not provide a good fit to SDSSJ2303+1422."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "residual_map = data - model_data\n",
+        "\n",
+        "plot_array(array=residual_map, title=\"Residual Map of fit to SDSSJ2303+1422\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Just like we did for the 1D `Gaussian` fitting examples, we can use the noise-map to compute the normalized residuals \n",
+        "and chi-squared map of the lens model."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "# The circular masking introduces zeros at the edge of the noise-map,\n",
+        "# which can lead to divide-by-zero errors.\n",
+        "# We set these values to 1.0e8, to ensure they do not contribute to the log likelihood.\n",
+        "noise_map_fit = noise_map\n",
+        "noise_map_fit[noise_map == 0] = 1.0e8\n",
+        "\n",
+        "normalized_residual_map = residual_map / noise_map_fit\n",
+        "\n",
+        "chi_squared_map = (normalized_residual_map) ** 2.0\n",
+        "\n",
+        "plot_array(\n",
+        "    array=normalized_residual_map,\n",
+        "    title=\"Normalized Residual Map of fit to SDSSJ2303+1422\",\n",
+        ")\n",
+        "plot_array(array=chi_squared_map, title=\"Chi Squared Map of fit to SDSSJ2303+1422\")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Finally, we can compute the `log_likelihood` of this lens model, which we will use in the next example to fit the \n",
+        "lens model to data with a non-linear search."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "chi_squared = np.sum(chi_squared_map)\n",
+        "noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))\n",
+        "\n",
+        "log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
+        "\n",
+        "print(log_likelihood)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Wrap Up__\n",
+        "\n",
+        "In this example, we introduced the astrophysical phenomena of strong gravitational lensing, and gave an overview of how\n",
+        "one can create a model for a strong lens system and fit it to imaging data. \n",
+        "\n",
+        "We ended by defining the log likelihood of the model-fit, which will form the `log_likelihood_function` of the\n",
+        "`Analysis` class we use in the next example, which fits this strong lens using **PyAutoFit**.\n",
+        "\n",
+        "There is one thing you should think about, how would we translate the above classes (e.g. `LightExponential`, \n",
+        "`MassIsothermal` and `Galaxy`) using the **PyAutoFit** `Model` and `Collection` objects? The `Galaxy` class contained \n",
+        "instances of the light and mass profile classes, meaning the standard use of the `Model` and `Collection` objects could \n",
+        "not handle this.\n",
+        "\n",
+        "This is where multi-level models come in, as will be shown in the next example!"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [],
+      "outputs": [],
+      "execution_count": null
+    }
+  ],
+  "metadata": {
+    "anaconda-cloud": {},
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.6.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/projects/cosmology/example_1_intro.py b/projects/cosmology/example_1_intro.py
new file mode 100644
index 0000000..f72f1de
--- /dev/null
+++ b/projects/cosmology/example_1_intro.py
@@ -0,0 +1,372 @@
+"""
+Project: Cosmology
+==================
+
+This project uses the astrophysical phenomena of Strong Gravitational Lensing to illustrate multi-level model
+composition and fitting with **PyAutoFit**.
+
+A strong gravitational lens is a system where two (or more) galaxies align perfectly down our line of sight from Earth
+such that the foreground galaxy's mass deflects the light of a background source galaxy(s).
+
+When the alignment is just right and the lens is massive enough, the background source galaxy appears multiple
+times. The schematic below shows such a system, where light-rays from the source are deflected around the lens galaxy
+to the observer following multiple distinct paths.
+
+![Schematic of Gravitational Lensing](https://raw.githubusercontent.com/Jammy2211/PyAutoLens/main/docs/overview/images/overview_1_lensing/schematic.jpg)
+**Credit: F. Courbin, S. G. Djorgovski, G. Meylan, et al., Caltech / EPFL / WMKO**
+https://www.cosmology.caltech.edu/~george/qsolens/
+
+As an observer, we don't see the source's true appearance (e.g. the red round blob of light). We only observe its
+light after it has been deflected and lensed by the foreground galaxies (e.g. as the two distinct red multiple images
+ in the image on the left). We also observe the emission of the foreground galaxy (in blue).
+
+You can read more about gravitational lensing as the following link:
+
+https://en.wikipedia.org/wiki/Gravitational_lens
+
+__PyAutoLens__
+
+Strong gravitational lensing is the original science case that sparked the development of **PyAutoFit**, which is
+a spin off of our astronomy software **PyAutoLens** `https://github.com/Jammy2211/PyAutoLens`.
+
+We'll use **PyAutoLens** to illustrate how the tools we developed with **PyAutoFit** allowed us to
+ensure **PyAutoLens**'s model fitting tools were extensible, easy to maintain and enabled intuitive model composition.
+
+__Multi-Level Models__
+
+Strong lensing is a great case study for using **PyAutoFit**, due to the multi-component nature of how one composes
+a strong lens model. A strong lens model consists of light and mass models of each galaxy in the lens system, where
+each galaxy is a model in itself. The galaxies are combined into one overall "lens model", which in later tutorials
+we will show may also have a Cosmological model.
+
+This example project uses **PyAutoFit** to compose and fit models of a strong lens, in particular highlighting
+**PyAutoFits** multi-level model composition.
+
+__Strong Lens Modeling__
+
+The models are fitted to Hubble Space Telescope imaging of a real strong lens system and will allow us to come up
+with a description of how light is deflected on its path through the Universe.
+
+This project consists of two example scripts / notebooks:
+
+ 1) `example_1_intro`: An introduction to strong lensing, and the various parts of the project's source code that are
+    used to represent a strong lens galaxy.
+
+ 2) `example_2_multi_level_model`: Using **PyAutoFit** to model a strong lens, with a strong emphasis on the
+    multi-level model API.
+
+__This Example__
+
+This introduction primarily focuses on what strong lensing is, how we define the individual model-components and fit
+a strong lens model to data. It does not make much use of **PyAutoFit**, but it does provide a clear understanding of
+the model so that **PyAutoFit**'s use in example 2 is clear.
+
+Note that import `import src as cosmo`. The package `src` contains all the code we need for this example Cosmology
+use case, and can be thought of as the source-code you would write to perform model-fitting via **PyAutoFit** for your
+problem of interest.
+
+The module `src/__init__.py` performs a series of imports that are used throughout this lecture to provide convenient
+access to different parts of the source code.
+"""
+
+# %matplotlib inline
+# from pyprojroot import here
+# workspace_path = str(here())
+# %cd $workspace_path
+# print(f"Working Directory has been set to `{workspace_path}`")
+
+import src as cosmo
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import signal
+from os import path
+
+"""
+__Plot__
+
+We will plot a lot of arrays of 2D data and grids of 2D coordinates in this example, so lets make a convenience 
+functions.
+"""
+
+
+def plot_array(array, title=None, norm=None):
+    plt.imshow(array, norm=norm)
+    plt.colorbar()
+    plt.title(title)
+    plt.show()
+    plt.close()
+
+
+def plot_grid(grid, title=None):
+    plt.scatter(x=grid[:, :, 0], y=grid[:, :, 1], s=1)
+    plt.title(title)
+    plt.show()
+    plt.close()
+
+
+"""
+__Data__
+
+First, lets load and plot Hubble Space Telescope imaging data of the strong gravitational lens called SDSSJ2303+1422, 
+where this data includes:
+ 
+ 1) The image of the strong lens, which is the data we'll fit.
+ 2) The noise in every pixel of this image, which will be used when evaluating the log likelihood.
+
+__Masking__
+
+When fitting 2D imaging data, it is common to apply a mask which removes regions of the image that are not relevant to
+the model fitting.
+
+For example, when fitting the strong lens, we remove the edges of the image where the lens and source galaxy's light is 
+not visible.
+
+In the strong lens image and noise map below, you can see this has already been performed, with the edge regions
+blank.
+"""
+dataset_path = path.join("projects", "cosmology", "dataset")
+
+data = np.load(file=path.join(dataset_path, "data.npy"))
+plot_array(array=data, title="Image of Strong Lens SDSSJ2303+1422")
+
+noise_map = np.load(file=path.join(dataset_path, "noise_map.npy"))
+plot_array(array=noise_map, title="Noise Map of Strong Lens SDSSJ2303+1422")
+
+"""
+In the image of the strong lens two distinct objects can be seen:
+
+ 1) A central blob of light, corresponding to the foreground lens galaxy whose mass is responsible for deflecting light.
+ 2) Two faint arcs of light in the bakcground, which is the lensed background source.
+ 
+__PSF__
+
+Another component of imaging data is the Point Spread Function (PSF), which describes how the light of the galaxies
+are blurred when they enter the Huble Space Telescope's. 
+
+This is because diffraction occurs when the light enters HST's optics, causing the light to smear out. The PSF is
+a two dimensional array that describes this blurring via a 2D convolution kernel.
+
+When fitting the data below and in the `log_likelihood_function`, you'll see that the PSF is used when creating the 
+model data. This is an example of how an `Analysis` class may be extended to include additional steps in the model
+fitting procedure.
+"""
+psf = np.load(file=path.join(dataset_path, "psf.npy"))
+plot_array(array=psf, title="Point Spread Function of Strong Lens SDSSJ2303+1422")
+
+
+""" 
+__Grid__
+ 
+To perform strong lensing, we need a grid of (x,y) coordinates which we map throughout the Universe as if their path
+is deflected. 
+
+For this, we create a simple 2D grid of coordinates below where the origin is (0.0, 0.0) and the size of
+a pixel is 0.05, which corresponds to the resolution of our image `data`. 
+
+This grid only contains (y,x) coordinates within the cricular mask that was applied to the data, as we only need to
+perform ray-tracing within this region.
+"""
+grid = np.load(file=path.join(dataset_path, "grid.npy"))
+
+plot_grid(
+    grid=grid,
+    title="Cartesian grid of (x,y) coordinates aligned with strong lens dataset",
+)
+
+"""
+__Light Profiles__
+
+Our model of a strong lens must include a description of the light of each galaxy, which we call a "light profile".
+In the source-code of this example project, specifically the module `src/light_profiles.py` you will see there
+are two light profile classes named `LightDeVaucouleurs` and `LightExponential`.
+
+These Python classes are the model components we will use to represent each galaxy's light and they behave analogous 
+to the `Gaussian` class seen in other tutorials. The input parameters of their `__init__` constructor (e.g. `centre`, 
+`axis_ratio`, `angle`) are their model parameters that may be fitted for by a non-linear search.
+
+These classes also contain functions which create an image from the light profile if an input grid of (x,y) 2D 
+coordinates are input, which we use below to create an image of a light profile.
+"""
+light_profile = cosmo.lp.LightExponential(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, intensity=1.0, effective_radius=2.0
+)
+light_image = light_profile.image_from_grid(grid=grid)
+
+plot_array(array=light_image, title="Image of an Exponential light profile.")
+
+"""
+__Mass Profiles__
+
+Our model also includes the mass of the foreground lens galaxy, called a 'mass profile'. In the source-code of the 
+example project, specifically the module `src/mass_profiles.py` you will see there is a mass profile class named 
+`MassIsothermal`. Like the light profile, this will be a model-component **PyAutoFit** fits via a non-linear search.
+
+The class also contains functions which create the "deflections angles", which describe the angles by which light is 
+deflected when it passes the mass of the foreground lens galaxy. These are subtracted from the (y,x) grid above to
+determine the original coordinates of the source galaxy before lensing.
+
+A higher mass galaxy, which bends light more, will have higher values of the deflection angles plotted below:
+"""
+mass_profile = cosmo.mp.MassIsothermal(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.5
+)
+mass_deflections = mass_profile.deflections_from_grid(grid=grid)
+
+plot_array(
+    array=mass_deflections[:, :, 0],
+    title="X-component of the deflection angles of a Isothermal mass profile.",
+)
+plot_array(
+    array=mass_deflections[:, :, 1],
+    title="Y-component of the deflection angles of a Isothermal mass profile.",
+)
+
+"""
+__Ray Tracing__
+
+The deflection angles describe how our (x,y) grid of coordinates are deflected by the mass of the foreground galaxy.
+
+We can therefore ray-trace the grid aligned with SDSSJ2303+1422 using the mass profile above and plot a grid of
+coordinates in the reference frame of before their light is gravitationally lensed:
+"""
+traced_grid = grid - mass_deflections
+
+plot_grid(grid=traced_grid, title="Cartesian grid of (x,y) traced coordinates.")
+
+"""
+By inputting this traced grid of (x,y) coordinates into our light profile, we can create an image of the galaxy as if
+it were gravitationally lensed by the mass profile.
+"""
+traced_light_image = light_profile.image_from_grid(grid=traced_grid)
+
+plot_array(
+    array=traced_light_image,
+    title="Image of a gravitationally lensed Exponential light profile.",
+)
+
+"""
+__Galaxy__
+
+In the `src/galaxy.py` module we define the `Galaxy` class, which is a collection of light and mass profiles at an 
+input redshift. For strong lens modeling, we have to use `Galaxy` objects, as the redshifts define how ray-tracing is
+performed.
+
+We create two instances of the `Galaxy` class, representing the lens and source galaxies in a strong lens system.
+"""
+light_profile = cosmo.lp.LightDeVaucouleurs(
+    centre=(0.01, 0.01), axis_ratio=0.9, angle=45.0, intensity=0.1, effective_radius=1.0
+)
+mass_profile = cosmo.mp.MassIsothermal(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.8
+)
+lens_galaxy = cosmo.Galaxy(
+    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]
+)
+
+light_profile = cosmo.lp.LightExponential(
+    centre=(0.1, 0.1), axis_ratio=0.5, angle=80.0, intensity=1.0, effective_radius=5.0
+)
+source_galaxy = cosmo.Galaxy(
+    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]
+)
+
+"""
+A galaxy's image is the sum of its light profile images, and its deflection angles are the sum of its mass profile
+deflection angles.
+
+To illustrate this, lets plot the lens galaxy's light profile image:
+"""
+galaxy_image = lens_galaxy.image_from_grid(grid=grid)
+
+plot_array(array=galaxy_image, title="Image of the Lens Galaxy.")
+
+"""
+__Data Fitting__
+
+We can create an overall image of the strong lens by:
+
+ 1) Creating an image of the lens galaxy.
+ 2) Computing the deflection angles of the lens galaxy.
+ 3) Ray-tracing light to the source galaxy reference frame and using its light profile to make its image.
+"""
+lens_image = lens_galaxy.image_from_grid(grid=grid)
+lens_deflections = lens_galaxy.deflections_from_grid(grid=grid)
+
+traced_grid = grid - lens_deflections
+
+source_image = source_galaxy.image_from_grid(grid=traced_grid)
+
+# The grid has zeros at its edges, which produce nans in the model image.
+# These lead to an ill-defined log likelihood, so we set them to zero.
+overall_image = np.nan_to_num(overall_image)
+
+plot_array(array=overall_image, title="Image of the overall Strong Lens System.")
+
+"""
+__Model Data__
+
+To produce the `model_data`, we now convolution the overall image with the Point Spread Function (PSF) of our
+observations. This blurs the image to simulate the telescope optics and pixelization used to observe the image.
+"""
+model_data = signal.convolve2d(overall_image, psf, mode="same")
+
+
+plot_array(array=model_data, title="Image of the overall Strong Lens System.")
+
+"""
+By subtracting this model image from the data, we can create a 2D residual map. This is equivalent to the residual maps
+we made in the 1D Gaussian examples, except for 2D imaging data.
+
+Clearly, the random lens model we used in this example does not provide a good fit to SDSSJ2303+1422.
+"""
+residual_map = data - model_data
+
+plot_array(array=residual_map, title="Residual Map of fit to SDSSJ2303+1422")
+
+"""
+Just like we did for the 1D `Gaussian` fitting examples, we can use the noise-map to compute the normalized residuals 
+and chi-squared map of the lens model.
+"""
+# The circular masking introduces zeros at the edge of the noise-map,
+# which can lead to divide-by-zero errors.
+# We set these values to 1.0e8, to ensure they do not contribute to the log likelihood.
+noise_map_fit = noise_map
+noise_map_fit[noise_map == 0] = 1.0e8
+
+normalized_residual_map = residual_map / noise_map_fit
+
+chi_squared_map = (normalized_residual_map) ** 2.0
+
+plot_array(
+    array=normalized_residual_map,
+    title="Normalized Residual Map of fit to SDSSJ2303+1422",
+)
+plot_array(array=chi_squared_map, title="Chi Squared Map of fit to SDSSJ2303+1422")
+
+"""
+Finally, we can compute the `log_likelihood` of this lens model, which we will use in the next example to fit the 
+lens model to data with a non-linear search.
+"""
+chi_squared = np.sum(chi_squared_map)
+noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))
+
+log_likelihood = -0.5 * (chi_squared + noise_normalization)
+
+print(log_likelihood)
+
+"""
+__Wrap Up__
+
+In this example, we introduced the astrophysical phenomena of strong gravitational lensing, and gave an overview of how
+one can create a model for a strong lens system and fit it to imaging data. 
+
+We ended by defining the log likelihood of the model-fit, which will form the `log_likelihood_function` of the
+`Analysis` class we use in the next example, which fits this strong lens using **PyAutoFit**.
+
+There is one thing you should think about, how would we translate the above classes (e.g. `LightExponential`, 
+`MassIsothermal` and `Galaxy`) using the **PyAutoFit** `Model` and `Collection` objects? The `Galaxy` class contained 
+instances of the light and mass profile classes, meaning the standard use of the `Model` and `Collection` objects could 
+not handle this.
+
+This is where multi-level models come in, as will be shown in the next example!
+"""
diff --git a/projects/cosmology/example_2_multi_level_model.ipynb b/projects/cosmology/example_2_multi_level_model.ipynb
new file mode 100644
index 0000000..39dc8e3
--- /dev/null
+++ b/projects/cosmology/example_2_multi_level_model.ipynb
@@ -0,0 +1,556 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Project: Cosmology\n",
+        "==================\n",
+        "\n",
+        "This project uses the astrophysical phenomena of Strong Gravitational Lensing to illustrate basic and advanced model\n",
+        "composition and fitting with **PyAutoFit**. The first tutorial described what a strong gravitational lens is and how\n",
+        "we build and fit a model of one.\n",
+        "\n",
+        "In this example, we use **PyAutoFit**'s multi-level models to compose a strong lens model consisting of a lens and\n",
+        "source galaxy, and fit it to the data on SDSSJ2303+1422.\n",
+        "\n",
+        "__Config Path__\n",
+        "\n",
+        "We first set up the path to this projects config files, which is located at `autofit_workspace/projects/cosmology/config`.\n",
+        "\n",
+        "This includes the default priors for the lens model, check it out!"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "import os\n",
+        "from os import path\n",
+        "from autoconf import conf\n",
+        "\n",
+        "cwd = os.getcwd()\n",
+        "config_path = path.join(cwd, \"projects\", \"cosmology\", \"config\")\n",
+        "conf.instance.push(new_path=config_path)\n",
+        "\n",
+        "# %matplotlib inline\n",
+        "# from pyprojroot import here\n",
+        "# workspace_path = str(here())\n",
+        "# %cd $workspace_path\n",
+        "# print(f\"Working Directory has been set to `{workspace_path}`\")\n",
+        "\n",
+        "import autofit as af\n",
+        "import src as cosmo\n",
+        "import matplotlib.pyplot as plt\n",
+        "import numpy as np"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Plot__\n",
+        "\n",
+        "First, lets again define the plotting convenience functions we used in the previous example."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "\n",
+        "def plot_array(array, title=None, norm=None):\n",
+        "    plt.imshow(array, norm=norm)\n",
+        "    plt.colorbar()\n",
+        "    plt.title(title)\n",
+        "    plt.show()\n",
+        "    plt.close()\n",
+        "\n",
+        "\n",
+        "def plot_grid(grid, title=None):\n",
+        "    plt.scatter(x=grid[:, :, 0], y=grid[:, :, 1], s=1)\n",
+        "    plt.title(title)\n",
+        "    plt.show()\n",
+        "    plt.close()\n"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Data__\n",
+        "\n",
+        "Now lets load and plot Hubble Space Telescope imaging data of the strong gravitational lens SDSSJ2303+1422."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "dataset_path = path.join(\"projects\", \"cosmology\", \"dataset\")\n",
+        "\n",
+        "data = np.load(file=path.join(dataset_path, \"data.npy\"))\n",
+        "plot_array(array=data, title=\"Image of Strong Lens SDSSJ2303+1422\")\n",
+        "\n",
+        "noise_map = np.load(file=path.join(dataset_path, \"noise_map.npy\"))\n",
+        "plot_array(array=noise_map, title=\"Noise Map of Strong Lens SDSSJ2303+1422\")\n",
+        "\n",
+        "psf = np.load(file=path.join(dataset_path, \"psf.npy\"))\n",
+        "plot_array(array=psf, title=\"Point Spread Function of Strong Lens SDSSJ2303+1422\")\n",
+        "\n",
+        "grid = np.load(file=path.join(dataset_path, \"grid.npy\"))\n",
+        "\n",
+        "plot_grid(\n",
+        "    grid=grid,\n",
+        "    title=\"Cartesian grid of (x,y) coordinates aligned with strong lens dataset\",\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Multi-level Model__\n",
+        "\n",
+        "In the previous example, we saw that we can use instances of the light profiles, mass profiles and galaxy objects to\n",
+        "perform strong lens ray-tracing calculations:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "light_profile = cosmo.lp.LightDeVaucouleurs(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, intensity=1.0, effective_radius=2.0\n",
+        ")\n",
+        "mass_profile = cosmo.mp.MassIsothermal(\n",
+        "    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.5\n",
+        ")\n",
+        "galaxy = cosmo.Galaxy(\n",
+        "    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]\n",
+        ")"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "In this example, we want to perform a model-fit using a non-linear search, where the `Galaxy` is a `Model`, but it\n",
+        "contains model subcomponents that are its individual light and mass profiles. \n",
+        "\n",
+        "Here is a pictoral representation of the model:\n",
+        "\n",
+        "![Strong Lens Model](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/lens_model.png?raw=true \"cluster\")\n",
+        "\n",
+        "__Model Composition__\n",
+        "\n",
+        "How do we compose a strong lens model where a `Galaxy` is a `Model`, but it contains the light and mass profiles \n",
+        "as `Model` themselves?\n",
+        "\n",
+        "We use **PyAutoFit**'s multi-level model composition:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "lens = af.Model(\n",
+        "    cls=cosmo.Galaxy,  # The overall model object uses this input.\n",
+        "    redshift=0.5,\n",
+        "    light_profile_list=[\n",
+        "        af.Model(cosmo.lp.LightDeVaucouleurs)\n",
+        "    ],  # These will be subcomponents of the model.\n",
+        "    mass_profile_list=[\n",
+        "        af.Model(cosmo.mp.MassIsothermal)\n",
+        "    ],  # These will be subcomponents of the model.\n",
+        ")\n",
+        "\n",
+        "print(lens.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Lets consider what is going on here:\n",
+        "\n",
+        " 1) We use a `Model` to create the overall model component. The `cls` input is the `Galaxy` class, therefore the \n",
+        "    overall model that is created is a `Galaxy`.\n",
+        "  \n",
+        " 2) **PyAutoFit** next inspects whether the key word argument inputs to the `Model` match any of the `__init__` \n",
+        "    constructor arguments of the `Galaxy` class. This determine if these inputs are to be composed as model \n",
+        "    subcomponents of the overall `Galaxy` model. \n",
+        "\n",
+        " 3) **PyAutoFit** matches the `light_profile_list` and  `mass_profile_list` inputs, noting they are passed as separate \n",
+        "    lists  containing the `LightDeVaucouleurs` and `MassIsothermal` class. They are both created as subcomponents of \n",
+        "    the overall `Galaxy` model.\n",
+        "  \n",
+        " 4) It also matches the `redshift` input, making it a fixed value of 0.5 for the model and not treating it as a \n",
+        "    free parameter.\n",
+        " \n",
+        "We can confirm this by printing the `total_free_parameters` of the lens, and noting it is 11 (6 parameters for \n",
+        "the `LightDeVaucouleurs` and 5 for the `MassIsothermal`)."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "print(lens.total_free_parameters)\n",
+        "print(lens.light_profile_list[0].total_free_parameters)\n",
+        "print(lens.mass_profile_list[0].total_free_parameters)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The `lens` behaves exactly like the model-components we are used to previously. For example, we can unpack its \n",
+        "individual parameters to customize the model, where below we:\n",
+        "\n",
+        " 1) Fix the light and mass profiles to the centre (0.0, 0.0).\n",
+        " 2) Customize the prior on the light profile `axis_ratio`.\n",
+        " 3) Fix the `axis_ratio` of the mass profile to 0.8."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "lens.light_profile_list[0].centre = (0.0, 0.0)\n",
+        "lens.light_profile_list[0].axis_ratio = af.UniformPrior(\n",
+        "    lower_limit=0.7, upper_limit=0.9\n",
+        ")\n",
+        "lens.light_profile_list[0].angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)\n",
+        "lens.light_profile_list[0].intensity = af.LogUniformPrior(\n",
+        "    lower_limit=1e-4, upper_limit=1e4\n",
+        ")\n",
+        "lens.light_profile_list[0].effective_radius = af.UniformPrior(\n",
+        "    lower_limit=0.0, upper_limit=5.0\n",
+        ")\n",
+        "\n",
+        "lens.mass_profile_list[0].centre = (0.0, 0.0)\n",
+        "lens.mass_profile_list[0].axis_ratio = 0.8\n",
+        "lens.mass_profile_list[0].angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)\n",
+        "lens.mass_profile_list[0].mass = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n",
+        "\n",
+        "print(lens.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Alternative API__\n",
+        "\n",
+        "We can create the `Galaxy` model component with the exact same customization by creating each profile as a `Model` and\n",
+        "passing these to the galaxy `Model`. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "light = af.Model(cosmo.lp.LightDeVaucouleurs)\n",
+        "\n",
+        "light.centre = af.UniformPrior(lower_limit=-0.05, upper_limit=0.05)\n",
+        "light.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=0.9)\n",
+        "light.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)\n",
+        "light.intensity = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4)\n",
+        "light.effective_radius = af.UniformPrior(lower_limit=0.0, upper_limit=5.0)\n",
+        "\n",
+        "\n",
+        "mass = af.Model(cosmo.mp.MassIsothermal)\n",
+        "\n",
+        "mass.centre = (0.0, 0.0)\n",
+        "mass.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=1.0)\n",
+        "mass.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)\n",
+        "mass.mass = af.UniformPrior(lower_limit=0.0, upper_limit=4.0)\n",
+        "\n",
+        "lens = af.Model(\n",
+        "    cosmo.Galaxy, redshift=0.5, light_profile_list=[light], mass_profile_list=[mass]\n",
+        ")\n",
+        "\n",
+        "print(lens.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can now create a model of our source galaxy using the same API."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "light = af.Model(cosmo.lp.LightExponential)\n",
+        "\n",
+        "light.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.3)\n",
+        "light.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.3)\n",
+        "light.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=1.0)\n",
+        "light.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)\n",
+        "light.intensity = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4)\n",
+        "light.effective_radius = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)\n",
+        "\n",
+        "source = af.Model(cosmo.Galaxy, redshift=1.0, light_profile_list=[light])\n",
+        "\n",
+        "print(source.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can now create our overall strong lens model, using a `Collection` in the same way we have seen previously."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "model = af.Collection(galaxies=af.Collection(lens=lens, source=source))\n",
+        "\n",
+        "print(model.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The model contains both galaxies in the strong lens, alongside all of their light and mass profiles.\n",
+        "\n",
+        "For every iteration of the non-linear search **PyAutoFit** generates an instance of this model, where all of the\n",
+        "`LightDeVaucouleurs`, `MassIsothermal` and `Galaxy` parameters of the are determined via their priors. \n",
+        "\n",
+        "An example instance is show below:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "instance = model.instance_from_prior_medians()\n",
+        "\n",
+        "print(\"Strong Lens Model Instance:\")\n",
+        "print(\"Lens Galaxy = \", instance.galaxies.lens)\n",
+        "print(\"Lens Galaxy Light = \", instance.galaxies.lens.profile_list)\n",
+        "print(\"Lens Galaxy Light Centre = \", instance.galaxies.lens.profile_list[0].centre)\n",
+        "print(\"Lens Galaxy Mass Centre = \", instance.galaxies.lens.mass_profile_list[0].centre)\n",
+        "print(\"Source Galaxy = \", instance.galaxies.source)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We have successfully composed a multi-level model, which we can fit via a non-linear search.\n",
+        "\n",
+        "At this point, you should check out the `Analysis` class of this example project, in the \n",
+        "module `projects/cosmology/src/analysis.py`. This class serves the same purpose that we have seen in the Gaussian 1D \n",
+        "examples, with the `log_likelihood_function` implementing the calculation we showed in the first tutorial.\n",
+        "\n",
+        "The `path_prefix1 and `name` inputs below sepciify the path and folder where the results of the model-fit are stored\n",
+        "in the output folder `autolens_workspace/output`. Results for this tutorial are writtent to hard-disk, due to the \n",
+        "longer run-times of the model-fit."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "search = af.DynestyStatic(\n",
+        "    path_prefix=path.join(\"projects\", \"cosmology\"),\n",
+        "    name=\"multi_level\",\n",
+        "    nlive=50,\n",
+        "    iterations_per_full_update=2500,\n",
+        ")\n",
+        "\n",
+        "analysis = cosmo.Analysis(data=data, noise_map=noise_map, psf=psf, grid=grid)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "If you comment out the code below, you will perform a lens model fit using the model and analysis class for \n",
+        "this project. However, this model-fit is slow to run, and it isn't paramount that you run it yourself.\n",
+        "\n",
+        "The animation below shows a slide-show of the lens modeling procedure. Many lens models are fitted to the data over\n",
+        "and over, gradually improving the quality of the fit to the data and looking more and more like the observed image.\n",
+        "\n",
+        "![Lens Modeling Animation](https://github.com/Jammy2211/auto_files/blob/main/lensmodel.gif?raw=true \"model\")"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "\n",
+        "result = search.fit(model=model, analysis=analysis)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "__Extensibility__\n",
+        "\n",
+        "This example project highlights how multi-level models can make certain model-fitting problem fully extensible. For \n",
+        "example:\n",
+        "\n",
+        " 1) A `Galaxy` class can be created using any combination of light and mass profiles, because it can wrap their\n",
+        " `image_from_grid` and `deflections_from_grid` methods as the sum of the individual profiles.\n",
+        " \n",
+        " 2) The overall strong lens model can contain any number of `Galaxy`'s, as their methods are used \n",
+        "    to implement the lensing calculations in the `Analysis` class and `log_likelihood_function`.\n",
+        "  \n",
+        "For problems of this nature, we can design and write code in a way that fully utilizes **PyAutoFit**'s multi-level\n",
+        "modeling features to compose and fits models of arbitrary complexity and dimensionality. \n",
+        "\n",
+        "__Galaxy Clusters__\n",
+        "\n",
+        "To illustrate this further, consider the following dataset which is called a \"strong lens galaxy cluster\":\n",
+        "\n",
+        "![Strong Lens Cluster](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/cluster_example.jpg?raw=true \"cluster\")\n",
+        "\n",
+        "For this strong lens, there are many tens of strong lens galaxies as well as multiple background source galaxies. \n",
+        "\n",
+        "However, despite it being a significantly more complex system than the single-galaxy strong lens we modeled above,\n",
+        "our use of multi-level models ensures that we can model such datasets without any additional code development, for\n",
+        "example:\n",
+        "\n",
+        "The lensing calculations in the source code `Analysis` object did not properly account for multiple galaxies \n",
+        "(called multi-plane ray tracing). This would need to be updated to properly model a galaxy cluster, but this\n",
+        "tutorial shows how a model can be composed for such a system."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "lens_0 = af.Model(\n",
+        "    cosmo.Galaxy,\n",
+        "    redshift=0.5,\n",
+        "    light_profile_list=[cosmo.lp.LightDeVaucouleurs],\n",
+        "    mass_profile_list=[cosmo.mp.MassIsothermal],\n",
+        ")\n",
+        "\n",
+        "lens_1 = af.Model(\n",
+        "    cosmo.Galaxy,\n",
+        "    redshift=0.5,\n",
+        "    light_profile_list=[cosmo.lp.LightDeVaucouleurs],\n",
+        "    mass_profile_list=[cosmo.mp.MassIsothermal],\n",
+        ")\n",
+        "\n",
+        "source_0 = af.Model(\n",
+        "    cosmo.Galaxy, redshift=1.0, light_profile_list=[af.Model(cosmo.lp.LightExponential)]\n",
+        ")\n",
+        "\n",
+        "# ... repeat for desired model complexity ...\n",
+        "\n",
+        "model = af.Collection(\n",
+        "    galaxies=af.Collection(\n",
+        "        lens_0=lens_0,\n",
+        "        lens_1=lens_1,\n",
+        "        source_0=source_0,\n",
+        "        # ... repeat for desired model complexity ...\n",
+        "    )\n",
+        ")\n",
+        "\n",
+        "print(model.info)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Here is a pictoral representation of a strong lens cluster as a multi-level model:\n",
+        "\n",
+        "![Strong Lens Cluster Model](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/lens_model_cluster.png?raw=true \"cluster\")\n",
+        "\n",
+        "__Wrap Up__\n",
+        "\n",
+        "Strong gravitational lensing is a great example of a problem that can be approached using multi-level models. \n",
+        "\n",
+        "At the core of this is how there are many different models one could imagine defining which describe the light or mass \n",
+        "of a galaxy. However, all of these models must derive the same fundamental property in order to fit the data, for\n",
+        "example the image of a light profile or the deflection angles of the mass profile.\n",
+        "\n",
+        "The multi-level nature of strong lensing is not unique, and is commonly found in my Astronomy problems and the \n",
+        "scientific literature in general. For example Astronomy problems:\n",
+        "\n",
+        " - Studies of galaxy structure, which represent the surface brightness distributions of galaxies as sums of Sersic\n",
+        " profiles (or other parametric equations) to quantify whether they are bulge-like or disk-like.\n",
+        " \n",
+        " - Studies of galaxy dynamics, which represent the mass distribution of galaxies as sums of profiles like the Isothermal\n",
+        " profile.\n",
+        " \n",
+        " - Studies of the activate galactic nuclei (AGN) of galaxies, where the different components of the AGN are represented\n",
+        " as different model components."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [],
+      "outputs": [],
+      "execution_count": null
+    }
+  ],
+  "metadata": {
+    "anaconda-cloud": {},
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.6.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/projects/cosmology/example_2_multi_level_model.py b/projects/cosmology/example_2_multi_level_model.py
new file mode 100644
index 0000000..e19f75f
--- /dev/null
+++ b/projects/cosmology/example_2_multi_level_model.py
@@ -0,0 +1,366 @@
+"""
+Project: Cosmology
+==================
+
+This project uses the astrophysical phenomena of Strong Gravitational Lensing to illustrate basic and advanced model
+composition and fitting with **PyAutoFit**. The first tutorial described what a strong gravitational lens is and how
+we build and fit a model of one.
+
+In this example, we use **PyAutoFit**'s multi-level models to compose a strong lens model consisting of a lens and
+source galaxy, and fit it to the data on SDSSJ2303+1422.
+
+__Config Path__
+
+We first set up the path to this projects config files, which is located at `autofit_workspace/projects/cosmology/config`.
+
+This includes the default priors for the lens model, check it out!
+"""
+
+import os
+from os import path
+from autoconf import conf
+
+cwd = os.getcwd()
+config_path = path.join(cwd, "projects", "cosmology", "config")
+conf.instance.push(new_path=config_path)
+
+# %matplotlib inline
+# from pyprojroot import here
+# workspace_path = str(here())
+# %cd $workspace_path
+# print(f"Working Directory has been set to `{workspace_path}`")
+
+import autofit as af
+import src as cosmo
+import matplotlib.pyplot as plt
+import numpy as np
+
+"""
+__Plot__
+
+First, lets again define the plotting convenience functions we used in the previous example.
+"""
+
+
+def plot_array(array, title=None, norm=None):
+    plt.imshow(array, norm=norm)
+    plt.colorbar()
+    plt.title(title)
+    plt.show()
+    plt.close()
+
+
+def plot_grid(grid, title=None):
+    plt.scatter(x=grid[:, :, 0], y=grid[:, :, 1], s=1)
+    plt.title(title)
+    plt.show()
+    plt.close()
+
+
+"""
+__Data__
+
+Now lets load and plot Hubble Space Telescope imaging data of the strong gravitational lens SDSSJ2303+1422.
+"""
+dataset_path = path.join("projects", "cosmology", "dataset")
+
+data = np.load(file=path.join(dataset_path, "data.npy"))
+plot_array(array=data, title="Image of Strong Lens SDSSJ2303+1422")
+
+noise_map = np.load(file=path.join(dataset_path, "noise_map.npy"))
+plot_array(array=noise_map, title="Noise Map of Strong Lens SDSSJ2303+1422")
+
+psf = np.load(file=path.join(dataset_path, "psf.npy"))
+plot_array(array=psf, title="Point Spread Function of Strong Lens SDSSJ2303+1422")
+
+grid = np.load(file=path.join(dataset_path, "grid.npy"))
+
+plot_grid(
+    grid=grid,
+    title="Cartesian grid of (x,y) coordinates aligned with strong lens dataset",
+)
+
+"""
+__Multi-level Model__
+
+In the previous example, we saw that we can use instances of the light profiles, mass profiles and galaxy objects to
+perform strong lens ray-tracing calculations:
+"""
+light_profile = cosmo.lp.LightDeVaucouleurs(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, intensity=1.0, effective_radius=2.0
+)
+mass_profile = cosmo.mp.MassIsothermal(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, mass=0.5
+)
+galaxy = cosmo.Galaxy(
+    redshift=0.5, light_profile_list=[light_profile], mass_profile_list=[mass_profile]
+)
+
+"""
+In this example, we want to perform a model-fit using a non-linear search, where the `Galaxy` is a `Model`, but it
+contains model subcomponents that are its individual light and mass profiles. 
+
+Here is a pictoral representation of the model:
+
+![Strong Lens Model](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/lens_model.png?raw=true "cluster")
+
+__Model Composition__
+
+How do we compose a strong lens model where a `Galaxy` is a `Model`, but it contains the light and mass profiles 
+as `Model` themselves?
+
+We use **PyAutoFit**'s multi-level model composition:
+"""
+lens = af.Model(
+    cls=cosmo.Galaxy,  # The overall model object uses this input.
+    redshift=0.5,
+    light_profile_list=[
+        af.Model(cosmo.lp.LightDeVaucouleurs)
+    ],  # These will be subcomponents of the model.
+    mass_profile_list=[
+        af.Model(cosmo.mp.MassIsothermal)
+    ],  # These will be subcomponents of the model.
+)
+
+print(lens.info)
+
+"""
+Lets consider what is going on here:
+
+ 1) We use a `Model` to create the overall model component. The `cls` input is the `Galaxy` class, therefore the 
+    overall model that is created is a `Galaxy`.
+  
+ 2) **PyAutoFit** next inspects whether the key word argument inputs to the `Model` match any of the `__init__` 
+    constructor arguments of the `Galaxy` class. This determine if these inputs are to be composed as model 
+    subcomponents of the overall `Galaxy` model. 
+
+ 3) **PyAutoFit** matches the `light_profile_list` and  `mass_profile_list` inputs, noting they are passed as separate 
+    lists  containing the `LightDeVaucouleurs` and `MassIsothermal` class. They are both created as subcomponents of 
+    the overall `Galaxy` model.
+  
+ 4) It also matches the `redshift` input, making it a fixed value of 0.5 for the model and not treating it as a 
+    free parameter.
+ 
+We can confirm this by printing the `total_free_parameters` of the lens, and noting it is 11 (6 parameters for 
+the `LightDeVaucouleurs` and 5 for the `MassIsothermal`).
+"""
+print(lens.total_free_parameters)
+print(lens.light_profile_list[0].total_free_parameters)
+print(lens.mass_profile_list[0].total_free_parameters)
+
+"""
+The `lens` behaves exactly like the model-components we are used to previously. For example, we can unpack its 
+individual parameters to customize the model, where below we:
+
+ 1) Fix the light and mass profiles to the centre (0.0, 0.0).
+ 2) Customize the prior on the light profile `axis_ratio`.
+ 3) Fix the `axis_ratio` of the mass profile to 0.8.
+"""
+
+lens.light_profile_list[0].centre = (0.0, 0.0)
+lens.light_profile_list[0].axis_ratio = af.UniformPrior(
+    lower_limit=0.7, upper_limit=0.9
+)
+lens.light_profile_list[0].angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)
+lens.light_profile_list[0].intensity = af.LogUniformPrior(
+    lower_limit=1e-4, upper_limit=1e4
+)
+lens.light_profile_list[0].effective_radius = af.UniformPrior(
+    lower_limit=0.0, upper_limit=5.0
+)
+
+lens.mass_profile_list[0].centre = (0.0, 0.0)
+lens.mass_profile_list[0].axis_ratio = 0.8
+lens.mass_profile_list[0].angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)
+lens.mass_profile_list[0].mass = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)
+
+print(lens.info)
+
+"""
+__Alternative API__
+
+We can create the `Galaxy` model component with the exact same customization by creating each profile as a `Model` and
+passing these to the galaxy `Model`. 
+"""
+light = af.Model(cosmo.lp.LightDeVaucouleurs)
+
+light.centre = af.UniformPrior(lower_limit=-0.05, upper_limit=0.05)
+light.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=0.9)
+light.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)
+light.intensity = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4)
+light.effective_radius = af.UniformPrior(lower_limit=0.0, upper_limit=5.0)
+
+
+mass = af.Model(cosmo.mp.MassIsothermal)
+
+mass.centre = (0.0, 0.0)
+mass.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=1.0)
+mass.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)
+mass.mass = af.UniformPrior(lower_limit=0.0, upper_limit=4.0)
+
+lens = af.Model(
+    cosmo.Galaxy, redshift=0.5, light_profile_list=[light], mass_profile_list=[mass]
+)
+
+print(lens.info)
+
+"""
+We can now create a model of our source galaxy using the same API.
+"""
+light = af.Model(cosmo.lp.LightExponential)
+
+light.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.3)
+light.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.3)
+light.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=1.0)
+light.angle = af.UniformPrior(lower_limit=0.0, upper_limit=180.0)
+light.intensity = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4)
+light.effective_radius = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)
+
+source = af.Model(cosmo.Galaxy, redshift=1.0, light_profile_list=[light])
+
+print(source.info)
+
+"""
+We can now create our overall strong lens model, using a `Collection` in the same way we have seen previously.
+"""
+model = af.Collection(galaxies=af.Collection(lens=lens, source=source))
+
+print(model.info)
+
+"""
+The model contains both galaxies in the strong lens, alongside all of their light and mass profiles.
+
+For every iteration of the non-linear search **PyAutoFit** generates an instance of this model, where all of the
+`LightDeVaucouleurs`, `MassIsothermal` and `Galaxy` parameters of the are determined via their priors. 
+
+An example instance is show below:
+"""
+instance = model.instance_from_prior_medians()
+
+print("Strong Lens Model Instance:")
+print("Lens Galaxy = ", instance.galaxies.lens)
+print("Lens Galaxy Light = ", instance.galaxies.lens.profile_list)
+print("Lens Galaxy Light Centre = ", instance.galaxies.lens.profile_list[0].centre)
+print("Lens Galaxy Mass Centre = ", instance.galaxies.lens.mass_profile_list[0].centre)
+print("Source Galaxy = ", instance.galaxies.source)
+
+"""
+We have successfully composed a multi-level model, which we can fit via a non-linear search.
+
+At this point, you should check out the `Analysis` class of this example project, in the 
+module `projects/cosmology/src/analysis.py`. This class serves the same purpose that we have seen in the Gaussian 1D 
+examples, with the `log_likelihood_function` implementing the calculation we showed in the first tutorial.
+
+The `path_prefix1 and `name` inputs below sepciify the path and folder where the results of the model-fit are stored
+in the output folder `autolens_workspace/output`. Results for this tutorial are writtent to hard-disk, due to the 
+longer run-times of the model-fit.
+"""
+
+search = af.DynestyStatic(
+    path_prefix=path.join("projects", "cosmology"),
+    name="multi_level",
+    nlive=50,
+    iterations_per_full_update=2500,
+)
+
+analysis = cosmo.Analysis(data=data, noise_map=noise_map, psf=psf, grid=grid)
+
+"""
+If you comment out the code below, you will perform a lens model fit using the model and analysis class for 
+this project. However, this model-fit is slow to run, and it isn't paramount that you run it yourself.
+
+The animation below shows a slide-show of the lens modeling procedure. Many lens models are fitted to the data over
+and over, gradually improving the quality of the fit to the data and looking more and more like the observed image.
+
+![Lens Modeling Animation](https://github.com/Jammy2211/auto_files/blob/main/lensmodel.gif?raw=true "model")
+"""
+
+result = search.fit(model=model, analysis=analysis)
+
+"""
+__Extensibility__
+
+This example project highlights how multi-level models can make certain model-fitting problem fully extensible. For 
+example:
+
+ 1) A `Galaxy` class can be created using any combination of light and mass profiles, because it can wrap their
+ `image_from_grid` and `deflections_from_grid` methods as the sum of the individual profiles.
+ 
+ 2) The overall strong lens model can contain any number of `Galaxy`'s, as their methods are used 
+    to implement the lensing calculations in the `Analysis` class and `log_likelihood_function`.
+  
+For problems of this nature, we can design and write code in a way that fully utilizes **PyAutoFit**'s multi-level
+modeling features to compose and fits models of arbitrary complexity and dimensionality. 
+
+__Galaxy Clusters__
+
+To illustrate this further, consider the following dataset which is called a "strong lens galaxy cluster":
+
+![Strong Lens Cluster](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/cluster_example.jpg?raw=true "cluster")
+
+For this strong lens, there are many tens of strong lens galaxies as well as multiple background source galaxies. 
+
+However, despite it being a significantly more complex system than the single-galaxy strong lens we modeled above,
+our use of multi-level models ensures that we can model such datasets without any additional code development, for
+example:
+
+The lensing calculations in the source code `Analysis` object did not properly account for multiple galaxies 
+(called multi-plane ray tracing). This would need to be updated to properly model a galaxy cluster, but this
+tutorial shows how a model can be composed for such a system.
+"""
+lens_0 = af.Model(
+    cosmo.Galaxy,
+    redshift=0.5,
+    light_profile_list=[cosmo.lp.LightDeVaucouleurs],
+    mass_profile_list=[cosmo.mp.MassIsothermal],
+)
+
+lens_1 = af.Model(
+    cosmo.Galaxy,
+    redshift=0.5,
+    light_profile_list=[cosmo.lp.LightDeVaucouleurs],
+    mass_profile_list=[cosmo.mp.MassIsothermal],
+)
+
+source_0 = af.Model(
+    cosmo.Galaxy, redshift=1.0, light_profile_list=[af.Model(cosmo.lp.LightExponential)]
+)
+
+# ... repeat for desired model complexity ...
+
+model = af.Collection(
+    galaxies=af.Collection(
+        lens_0=lens_0,
+        lens_1=lens_1,
+        source_0=source_0,
+        # ... repeat for desired model complexity ...
+    )
+)
+
+print(model.info)
+
+"""
+Here is a pictoral representation of a strong lens cluster as a multi-level model:
+
+![Strong Lens Cluster Model](https://github.com/rhayes777/PyAutoFit/blob/main/docs/overview/image/lens_model_cluster.png?raw=true "cluster")
+
+__Wrap Up__
+
+Strong gravitational lensing is a great example of a problem that can be approached using multi-level models. 
+
+At the core of this is how there are many different models one could imagine defining which describe the light or mass 
+of a galaxy. However, all of these models must derive the same fundamental property in order to fit the data, for
+example the image of a light profile or the deflection angles of the mass profile.
+
+The multi-level nature of strong lensing is not unique, and is commonly found in my Astronomy problems and the 
+scientific literature in general. For example Astronomy problems:
+
+ - Studies of galaxy structure, which represent the surface brightness distributions of galaxies as sums of Sersic
+ profiles (or other parametric equations) to quantify whether they are bulge-like or disk-like.
+ 
+ - Studies of galaxy dynamics, which represent the mass distribution of galaxies as sums of profiles like the Isothermal
+ profile.
+ 
+ - Studies of the activate galactic nuclei (AGN) of galaxies, where the different components of the AGN are represented
+ as different model components.
+"""
diff --git a/projects/cosmology/src/__init__.py b/projects/cosmology/src/__init__.py
new file mode 100644
index 0000000..217c871
--- /dev/null
+++ b/projects/cosmology/src/__init__.py
@@ -0,0 +1,4 @@
+from src import light_profiles as lp
+from src import mass_profiles as mp
+from src.galaxy import Galaxy
+from src.analysis import Analysis
diff --git a/projects/cosmology/src/analysis.py b/projects/cosmology/src/analysis.py
new file mode 100644
index 0000000..6a4b16e
--- /dev/null
+++ b/projects/cosmology/src/analysis.py
@@ -0,0 +1,248 @@
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import signal
+from typing import List
+
+import autofit as af
+
+
+class Analysis(af.Analysis):
+    def __init__(
+        self, data: np.ndarray, noise_map: np.ndarray, psf: np.ndarray, grid: np.ndarray
+    ):
+        """
+        The analysis class for the **PyAutoFit** example Astronomy project on gravitational lensing.
+
+        This class contains imaging data of a strong gravitational lens, and it fits it with a lens model which can
+        include:
+
+        1) The light and mass of the foreground galaxy responsible for strong lensing.
+        2) The light of the background source galaxy which is observed multiple times.
+        3) The light and mass of additional galaxies nearby, that may be included in the strong lens model.
+
+        The imaging data contains:
+
+        1) An image of the strong lens.
+        2) The noise in every pixel of that image.
+        3) The Point Spread Function (PSF) describing how the optics of the telescope blur the image.
+        4) A (y,x) grid of coordinates describing the locations of the image pixels in a unit system, which are used
+           to ray-trace the strong lens model.
+
+        This project is a scaled down version of the Astronomy project **PyAutoLens**, which was the original project
+        from which PyAutoFit is an offshoot!
+
+        https://github.com/Jammy2211/PyAutoLens
+
+        Parameters
+        ----------
+        data
+            The image containing the observation of the strong gravitational lens that is fitted.
+        noise_map
+            The RMS noise values of the image data, which is folded into the log likelihood calculation.
+        grid
+            The (y, x) coordinates of the image from which the lensing calculation is performed and model image is
+            computed using.
+        """
+
+        self.data = data
+        self.noise_map = noise_map
+        self.psf = psf
+        self.grid = grid
+
+        # The circular masking introduces zeros at the edge of the noise-map,
+        # which can lead to divide-by-zero errors.
+        # We set these values to 1.0e8, to ensure they do not contribute to the log likelihood.
+        self.noise_map_fit = noise_map
+        self.noise_map_fit[noise_map == 0] = 1.0e8
+
+    def log_likelihood_function(self, instance) -> float:
+        """
+        The `log_likelihood_function` of the strong lensing example, which performs the following step:
+
+        1) Using the lens model passed into the function (whose parameters are set via the non-linear search) create
+        an image of this strong lens. This uses gravitational lensing ray-tracing calculations.
+
+        2) Subtract this model image from the data and compute its residuals, chi-squared and likelihood via the
+        noise map.
+
+        Parameters
+        ----------
+        instance
+            An instance of the lens model set via the non-linear search.
+
+        Returns
+        -------
+        float
+            The log likelihood value of this particular lens model.
+        """
+
+        """
+        The 'instance' that comes into this method contains the `Galaxy`'s we setup in the `Model` and `Collection`,
+        which can be seen by uncommenting the code below.
+        """
+
+        # print("Lens Model Instance:")
+        # print("Lens Galaxy = ", instance.galaxies.lens)
+        # print("Lens Galaxy Bulge = ", instance.galaxies.lens.light_profile_list)
+        # print("Lens Galaxy Bulge Centre = ", instance.galaxies.lens.light_profile_list[0].centre)
+        # print("Lens Galaxy Mass Centre = ", instance.galaxies.lens.mass_profile_list[0].centre)
+        # print("Source Galaxy = ", instance.galaxies.source)
+
+        """
+        The methods of the `Galaxy` class are available, making it easy to fit
+        the lens model.
+        """
+
+        model_data = self.model_data_from_instance(instance=instance)
+
+        residual_map = self.data - model_data
+
+        normalized_residual_map = residual_map / self.noise_map_fit
+
+        chi_squared_map = (normalized_residual_map) ** 2.0
+
+        chi_squared = np.sum(chi_squared_map)
+        noise_normalization = np.sum(np.log(2 * np.pi * self.noise_map**2.0))
+
+        log_likelihood = -0.5 * (chi_squared + noise_normalization)
+
+        return log_likelihood
+
+    def visualize(self, paths: af.DirectoryPaths, instance, during_analysis):
+        """
+        Visualizes the maximum log likelihood model-fit to the strong lens data so far in the non-linear search, as
+        well as at the end of the search.
+
+        Parameters
+        ----------
+        paths
+            The **PyAutoFit** `Paths` object which includes the path to the output image folder.
+        instance
+            An instance of the lens model set via the non-linear search.
+        during_analysis
+            If the visualization is being performed during the non-linear search or one it is complete.
+        """
+
+        model_data = self.model_data_from_instance(instance=instance)
+
+        residual_map = self.data - model_data
+        normalized_residual_map = residual_map / self.noise_map_fit
+        chi_squared_map = (normalized_residual_map) ** 2.0
+
+        def plot_array(array, name, title=None, norm=None):
+            import os
+
+            os.makedirs(paths.image_path, exist_ok=True)
+
+            plt.imshow(array, norm=norm)
+            plt.colorbar()
+            plt.title(title)
+            plt.savefig(f"{paths.image_path}/{name}.png")
+            plt.close()
+
+        plot_array(array=self.data, title="Data", name="data")
+        plot_array(array=self.noise_map, title="Noise Map", name="noise_map")
+        plot_array(array=model_data, title="Model Data", name="model_data")
+        plot_array(array=residual_map, title="Residual Map", name="residual_map")
+        plot_array(
+            array=normalized_residual_map,
+            title="Normalized Residual Map",
+            name="normalized_residual_map",
+        )
+        plot_array(
+            array=chi_squared_map, title="Chi-Squared Map", name="chi_squared_map"
+        )
+
+    def traced_grid_2d_from(self, instance) -> List[np.ndarray]:
+        """
+        This function performs ray-tracing calculations describing how light is gravitationally lensed
+        on its path through the Universe by a foreground galaxy.
+
+        For the purpose of illustrating **PyAutoFit** you do not need to understand what this function is
+        doing, the main thing to note is that it allows us to create a model lensed image of a strong
+        lens model.
+
+        Nevertheless, if you are curious, a simple description of the steps performed is as follows:
+
+        1) Go to the foreground `lens` galaxy in the model `instance` (the lower redshift galaxy) and use
+           its mass profiles to compute the deflection angles of light-rays, describing how the light
+           of the galaxies behind it are deflected.
+
+        2) Subtract these deflection angles from the image-plane coordinates to determine where these
+           light-rays end up after bending, meaning that they tell us how the deflected image of the
+           background galaxies appear.
+
+        3) Repeat this galaxy by galaxy, working out how they deflect light and how this changes the
+           appearance of the images of the background galaxies.
+
+        There are strong lens systems in nature which consist of multiple galaxies at different redshifts
+        (e.g. distances from Earth). Such systems require their own bespoke ray-tracing calculations,
+        called "multi-plane" ray-tracing. This function does not perform these calculations, as they are
+        quite complicated and would make the code below less clear.
+
+        Parameters
+        ----------
+        instance
+            An instance of the lens model set via the non-linear search.
+
+        Returns
+        -------
+        grids_at_planes
+            The deflected (y,x) coordinates of the original image at every plane of the strong lens system, which can
+            be used for computing its image.
+        """
+
+        lens = instance.galaxies[0]
+
+        deflections = lens.deflections_from_grid(grid=self.grid)
+
+        lensed_grid = self.grid - deflections
+
+        return lensed_grid
+
+    def model_data_from_instance(self, instance):
+        """
+        Create the image of a strong lens system, accounting for the complicated ray-tracing calculations that
+        describe how light is gravitationally lensed on its path through the Universe by nearby galaxies.
+
+        For the purpose of illustrating **PyAutoFit** you do not need to understand what this function is doing, the
+        main thing to note is that it allows us to create a model lensed image of a strong lens model.
+
+        Nevertheless, if you are curious, it simply computes deflected (y,x) Cartesian grids by performing ray-tracing
+        using the `grids_via_tracer_from` function and uses each galaxy's light profile, based on where these
+        deflected grids fall on the image, to evaluate the appearance of each galaxy. This is summed to create the
+        overall image of the strong lens system.
+
+        NOTE: To avoid an AstroPy dependency the calculation below is not strictly correct, as it omits a Cosmology
+        calculation that rescales the deflection of light. This does not impact the project's
+        illustration of **PyAutoFit**.
+
+        Parameters
+        ----------
+        instance
+            An instance of the lens model set via the non-linear search.
+
+        Returns
+        -------
+        image
+            The image of this strong lens model.
+        """
+
+        lens = instance.galaxies[0]
+        lens_image = lens.image_from_grid(grid=self.grid)
+
+        source = instance.galaxies[1]
+        traced_grid = self.traced_grid_2d_from(instance=instance)
+        source_image = source.image_from_grid(grid=traced_grid)
+
+        overall_image = lens_image + source_image
+
+        # The grid has zeros at its edges, which produce nans in the model image.
+        # These lead to an ill-defined log likelihood, so we set them to zero.
+        overall_image = np.nan_to_num(overall_image)
+
+        model_data = signal.convolve2d(overall_image, self.psf, mode="same")
+
+        return model_data
diff --git a/projects/cosmology/src/galaxy.py b/projects/cosmology/src/galaxy.py
new file mode 100644
index 0000000..9fd5aa0
--- /dev/null
+++ b/projects/cosmology/src/galaxy.py
@@ -0,0 +1,74 @@
+import numpy as np
+from typing import Optional, List
+
+
+class Galaxy:
+    def __init__(
+        self,
+        redshift: float,
+        light_profile_list: Optional[List] = None,
+        mass_profile_list: Optional[List] = None,
+    ):
+        """
+        A galaxy, which contains light and mass profiles at a specified redshift.
+
+        Parameters
+        ----------
+        redshift
+            The redshift of the galaxy.
+        light_profile_list
+            A list of the galaxy's light profiles.
+        mass_profile_list
+            A list of the galaxy's mass profiles.
+        """
+
+        self.redshift = redshift
+        self.light_profile_list = light_profile_list
+        self.mass_profile_list = mass_profile_list
+
+    def image_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the summed 2D image of all of the galaxy's light profiles using an input
+        grid of Cartesian (y,x) coordinates.
+
+        If the galaxy has no light profiles, a grid of zeros is returned.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the original reference frame of the grid.
+        """
+        if self.light_profile_list is not None:
+            return sum(
+                map(lambda p: p.image_from_grid(grid=grid), self.light_profile_list)
+            )
+        return np.zeros((grid.shape[0],))
+
+    def deflections_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the summed (y,x) deflection angles of the galaxy's mass profiles
+        using a grid of Cartesian (y,x) coordinates.
+
+        If the galaxy has no mass profiles, two grid of zeros are returned.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the original reference frame of the grid.
+        """
+        if self.mass_profile_list is not None:
+            return sum(
+                map(
+                    lambda p: p.deflections_from_grid(grid=grid), self.mass_profile_list
+                )
+            )
+        return np.zeros((grid.shape[0], 2))
+
+
+class Redshift(float):
+    def __new__(cls, redshift):
+        # noinspection PyArgumentList
+        return float.__new__(cls, redshift)
+
+    def __init__(self, redshift):
+        float.__init__(redshift)
diff --git a/projects/cosmology/src/geometry_profiles.py b/projects/cosmology/src/geometry_profiles.py
new file mode 100644
index 0000000..674fa96
--- /dev/null
+++ b/projects/cosmology/src/geometry_profiles.py
@@ -0,0 +1,118 @@
+import typing
+import numpy as np
+
+
+class GeometryProfile:
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+    ):
+        """
+        Abstract base class for the geometry of a profile representing the light or mass of a galaxy.
+
+        Using the centre, axis-ratio and position angle of the profile this class describes how to convert a
+        (y,x) grid of Cartesian coordinates to the elliptical geometry of the profile.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        """
+
+        self.centre = centre
+        self.axis_ratio = axis_ratio
+        self.angle = angle
+
+    def transformed_to_reference_frame_grid_from(self, grid: np.ndarray):
+        """
+        Transform a grid of (y,x) coordinates to the geometric reference frame of the profile via a translation using
+        its `centre` and a rotation using its `angle`.
+
+        This performs the following steps:
+
+        1) Translate the input (y,x) coordinates from the (0.0, 0.0) origin to the centre of the profile by subtracting
+           the profile centre.
+
+        2) Compute the radial distance of every translated coordinate from the centre.
+
+        3) Rotate the coordinates from the above step counter-clockwise from the positive x-axis by the
+           profile's `angle`.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinate grid in its original reference frame.
+        """
+
+        shifted_grid = grid - self.centre
+        effective_radius = (
+            (shifted_grid[:, :, 0] ** 2.0 + shifted_grid[:, :, 1] ** 2.0)
+        ) ** 0.5
+
+        theta_coordinate_to_profile = np.arctan2(
+            shifted_grid[:, :, 1], shifted_grid[:, :, 0]
+        ) - np.radians(self.angle)
+
+        transformed_grid = np.zeros(grid.shape)
+
+        transformed_grid[:, :, 0] = effective_radius * np.cos(
+            theta_coordinate_to_profile
+        )
+        transformed_grid[:, :, 1] = effective_radius * np.sin(
+            theta_coordinate_to_profile
+        )
+
+        return transformed_grid
+
+    def rotated_grid_from_reference_frame_from(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Rotate a grid of (y,x) coordinates which have been transformed to the elliptical reference frame of a profile
+        back to the original unrotated coordinate frame.
+
+        This performs the following steps:
+
+        1) Rotate the coordinates from the elliptical reference frame back to the original unrotated coordinate frame
+           using the profile's `angle`.
+
+        2) Translate the coordinates from the centre of the profile back to the (0.0, 0.0) origin by adding the
+        profile centre.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the reference frame of an elliptical profile.
+        """
+        cos_angle = np.cos(np.radians(self.angle))
+        sin_angle = np.sin(np.radians(self.angle))
+
+        transformed_grid = np.zeros(grid.shape)
+
+        transformed_grid[:, :, 0] = np.add(
+            np.multiply(grid[:, :, 0], cos_angle),
+            -np.multiply(grid[:, :, 1], sin_angle),
+        )
+        transformed_grid[:, :, 1] = np.add(
+            np.multiply(grid[:, :, 0], sin_angle), np.multiply(grid[:, :, 1], cos_angle)
+        )
+
+        return transformed_grid
+
+    def elliptical_radii_grid_from(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Convert a grid of (y,x) coordinates to a grid of elliptical radii.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the reference frame of the elliptical profile.
+        """
+
+        return (
+            (grid[:, :, 1] ** 2.0) + (grid[:, :, 0] / self.axis_ratio) ** 2.0
+        ) ** 0.5
diff --git a/projects/cosmology/src/light_profiles.py b/projects/cosmology/src/light_profiles.py
new file mode 100644
index 0000000..dc5ae58
--- /dev/null
+++ b/projects/cosmology/src/light_profiles.py
@@ -0,0 +1,142 @@
+from src import geometry_profiles as gp
+
+import typing
+import numpy as np
+
+
+class LightProfile(gp.GeometryProfile):
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        Abstract base class for a light profile, which describes the emission of a galaxy as a
+        function of radius.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(centre=centre, axis_ratio=axis_ratio, angle=angle)
+
+        self.intensity = intensity
+        self.effective_radius = effective_radius
+
+
+class LightDeVaucouleurs(LightProfile):
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        The De Vaucouleurs light profile often used in Astronomy to represent the bulge of galaxies.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(
+            centre=centre,
+            axis_ratio=axis_ratio,
+            angle=angle,
+            intensity=intensity,
+            effective_radius=effective_radius,
+        )
+
+    def image_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the image of the De Vaucouleurs light profile on a grid of Cartesian (y,x) coordinates, which are
+        first translated to the profile's reference frame.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates where the image is computed.
+        """
+        grid_transformed = self.transformed_to_reference_frame_grid_from(grid=grid)
+        grid_elliptical_radii = self.elliptical_radii_grid_from(grid=grid_transformed)
+
+        return self.intensity * np.exp(
+            -7.66924
+            * ((grid_elliptical_radii / self.effective_radius) ** (1.0 / 7.66924) - 1.0)
+        )
+
+
+class LightExponential(LightProfile):
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        The Exponential light profile often used in Astronomy to represent the disk of galaxies.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(
+            centre=centre,
+            axis_ratio=axis_ratio,
+            angle=angle,
+            intensity=intensity,
+            effective_radius=effective_radius,
+        )
+
+    def image_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the image of the light profile on a grid of Cartesian (y,x) coordinates.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates where the image is computed.
+        """
+        grid_transformed = self.transformed_to_reference_frame_grid_from(grid=grid)
+        grid_elliptical_radii = self.elliptical_radii_grid_from(grid=grid_transformed)
+
+        return self.intensity * np.exp(
+            -1.67838
+            * ((grid_elliptical_radii / self.effective_radius) ** (1.0 / 1.67838) - 1.0)
+        )
diff --git a/projects/cosmology/src/mass_profiles.py b/projects/cosmology/src/mass_profiles.py
new file mode 100644
index 0000000..0bfa387
--- /dev/null
+++ b/projects/cosmology/src/mass_profiles.py
@@ -0,0 +1,109 @@
+from src import geometry_profiles as gp
+
+import typing
+import numpy as np
+
+
+class MassProfile(gp.GeometryProfile):
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        mass: float = 1.0,
+    ):
+        """
+        Abstract base class for a mass profile, which describes the mass of a galaxy as a function of radius.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        mass
+            The mass intensity of the profile, which is the Einstein effective_radius in arc-seconds.
+        """
+
+        super().__init__(centre=centre, axis_ratio=axis_ratio, angle=angle)
+
+        self.mass = mass
+
+
+class MassIsothermal(MassProfile):
+    def __init__(
+        self,
+        centre: typing.Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        mass: float = 1.0,
+    ):
+        """
+        The elliptical isothermal mass distribution often used in Astronomy to represent the combined mass of stars
+        and dark matter in galaxies.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        mass
+            The mass intensity of the profile, which is the Einstein effective_radius in arc-seconds.
+        """
+        super().__init__(centre=centre, axis_ratio=axis_ratio, angle=angle, mass=mass)
+
+    def psi_from(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns `psi`, a value required when computing the deflections of the elliptical isothermal mass distribution,
+        where:
+
+        psi = sqrt((q^2 * x^2) + y^2)
+
+        Parameters
+        ----------
+        grid
+            The (y,x) coordinates of the grid where the `psi`'s are computed.
+        """
+
+        return (
+            (self.axis_ratio**2.0 * grid[:, :, 0] ** 2.0) + grid[:, :, 1] ** 2.0
+        ) ** 0.5
+
+    def deflections_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the deflection angles on a grid of (y,x) coordinates, which describe how
+        the isothermal mass profile bends light via the effect of gravitational lensing.
+
+        Parameters
+        ----------
+        grid
+            The grid of (y,x) arc-second coordinates the deflection angles are computed on.
+        """
+
+        grid_transformed = self.transformed_to_reference_frame_grid_from(grid=grid)
+
+        factor = 2.0 * self.mass * self.axis_ratio / np.sqrt(1 - self.axis_ratio**2)
+
+        psi = self.psi_from(grid=grid_transformed)
+
+        deflections = np.zeros(grid.shape)
+
+        deflections[:, :, 0] = factor * np.arctan(
+            np.divide(
+                np.multiply(np.sqrt(1 - self.axis_ratio**2), grid_transformed[:, :, 0]),
+                psi,
+            )
+        )
+        deflections[:, :, 1] = factor * np.arctanh(
+            np.divide(
+                np.multiply(np.sqrt(1 - self.axis_ratio**2), grid_transformed[:, :, 1]),
+                psi,
+            )
+        )
+
+        return self.rotated_grid_from_reference_frame_from(grid=deflections)
diff --git a/scripts/howtofit/chapter_1_introduction/tutorial_8_astronomy_example.py b/scripts/howtofit/chapter_1_introduction/tutorial_8_astronomy_example.py
new file mode 100644
index 0000000..bb50b7c
--- /dev/null
+++ b/scripts/howtofit/chapter_1_introduction/tutorial_8_astronomy_example.py
@@ -0,0 +1,1063 @@
+"""
+Tutorial 8: Astronomy Example
+=============================
+
+In this tutorial, we'll apply the tools we've learned in this chapter to tackle a read world problem and fit
+2D Hubble Space Telescope imaging data of galaxies.
+
+One of the most well-known results in astronomy is the Hubble sequence, which classifies galaxies based on their
+visual appearance. The sequence is divided into three broad groups:
+
+- **Early-type**: Galaxies that are round in shape with a smooth blob of light, often called a "bulge".
+
+- **Late-type**: Galaxies that are elliptical in shape with a flattened light distribution, often called a "disk".
+
+- **Irregular**: Galaxies that do not have a regular shape.
+
+Here is the Hubble Sequence, showing early-type galaxies on the left and late-type galaxies on the right:
+
+![Hubble Sequence](https://github.com/Jammy2211/autofit_workspace/blob/main/scripts/howtofit/chapter_1_introduction/images/hubble_sequence.png?raw=true)
+
+This example fits a Hubble Space Telescope image of a galaxy with two models: one representing bulge-like early-type
+galaxies and one representing disk-like late-type galaxies. This will help us determine whether the galaxy is an early
+or late-type galaxy.
+
+The aim of this example is to illustrate that everything you have learnt in this chapter can be applied to real-world
+problems, and to show you how many of the tools and challenges you have learnt about are used in practice.
+
+After fitting the Hubble Space Telescope data, we will explore some common model-fitting problems astronomers face
+when classifying galaxies and measuring their properties. This will be an open exercise for you to consider how these
+problems can be addressed.
+
+__Overview__
+
+In this tutorial, we will:
+
+- Use the tools from tutorials 1, 2, and 3 to set up a model that fits 2D astronomical images of galaxies with light
+  profile models.
+
+- Encounter and discuss the challenges of model fitting described in tutorial 4, along with strategies to overcome
+  these challenges.
+
+- Use the result properties introduced in tutorial 5 to interpret the model fit and determine whether the galaxy is an
+  early or late-type galaxy.
+
+ __Contents__
+
+- **Plot**: Convensional plotting functions for 2D data and grids of 2D coordinates.
+- **Data**: Load and plot Hubble Space Telescope imaging data of a galaxy.
+- **Mask**: Apply a mask to the data to remove regions of the image that are not relevant to the model fitting.
+- **PSF**: Load and plot the Point Spread Function (PSF) of the telescope.
+- **Grid**: Create a grid of (y,x) coordinates that overlap the observed galaxy data.
+- **Light Profiles**: Define light profile classes representing the light of galaxies.
+- **Model Data**: Create the model image of a galaxy by convolving its light profile with the PSF.
+- **Model**: Define a model with a light profile to fit the galaxy data.
+- **Analysis**: Define the log likelihood function, which compares the model image to the observed image.
+- **Model Fit**: Fit the model to the data and display the results.
+- **Result**: Interpret the model fit to determine whether the galaxy is an early or late-type galaxy.
+- **Bulgey**: Repeat the fit using a bulgey light profile to determine the galaxy's type.
+- **Model Mismatch**: Analyze the challenges from model mismatches in galaxy classification.
+- **Extensions**: Illustrate examples of how this problem can be extended and the challenges that arise.
+- **Chapter Wrap Up**: Summarize the completion of Chapter 1 and its applications to real astronomy.
+"""
+
+# from autoconf import setup_notebook; setup_notebook()
+
+from os import path
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import signal
+from typing import Tuple
+
+import autofit as af
+
+"""
+__Plot__
+
+We will plot a lot of arrays of 2D data and grids of 2D coordinates in this example, so lets make a convenience 
+functions.
+"""
+
+
+def plot_array(array, title=None, norm=None, filename=None):
+    plt.imshow(array, norm=norm)
+    plt.colorbar()
+    plt.title(title)
+    if filename is not None:
+        plt.savefig(filename)
+    plt.show()
+    plt.clf()
+    plt.close()
+
+
+def plot_grid(grid, title=None):
+    plt.scatter(x=grid[:, :, 0], y=grid[:, :, 1], s=1)
+    plt.title(title)
+    plt.show()
+    plt.clf()
+    plt.close()
+
+
+"""
+__Data__
+
+First, let's load and plot Hubble Space Telescope imaging data of a galaxy. This data includes:
+
+1) The image of the galaxy, which is the data we'll fit.
+2) The noise in each pixel of this image, which will be used to evaluate the log likelihood.
+
+The noise-map has a few strange off-centre features which are an artefact of the telescope. Don't worry about these
+features.
+"""
+dataset_path = path.join("dataset", "howtofit", "chapter_1", "astro", "simple")
+
+"""
+__Dataset Auto-Simulation__
+
+If the dataset does not already exist on your system, it will be created by running the corresponding
+simulator script. This ensures that all example scripts can be run without manually simulating data first.
+"""
+if not path.exists(dataset_path):
+    import subprocess
+    import sys
+
+    subprocess.run(
+        [sys.executable, "scripts/simulators/simulators.py"],
+        check=True,
+    )
+
+data = np.load(file=path.join(dataset_path, "data.npy"))
+plot_array(array=data, title="Image of Galaxy")
+
+noise_map = np.load(file=path.join(dataset_path, "noise_map.npy"))
+plot_array(array=noise_map, title="Noise Map of Galaxy")
+
+"""
+__Mask__
+
+When fitting 2D imaging data, it is common to apply a mask which removes regions of the image that are not relevant to
+the model fitting.
+
+For example, when fitting the galaxy, we remove the edges of the image where the galaxy's light is not visible.
+
+We load and plot the mask below to show you how it is applied to the data, and we will use it in 
+the `log_likelihood_function` below to ensure these regions are not fitted.
+"""
+mask = np.load(file=path.join(dataset_path, "mask.npy"))
+plot_array(array=mask, title="Mask of Galaxy")
+
+"""
+In the image of the galaxy a bright blob of light is clearly visible, which is the galaxy we'll fit with a model.
+
+__PSF__
+
+Another important component of imaging data is the Point Spread Function (PSF), which describes how the light from 
+galaxies is blurred when it enters the Hubble Space Telescope.
+
+This blurring occurs due to diffraction as light passes through the HST's optics. The PSF is represented as a 
+two-dimensional array, which acts as a 2D convolution kernel.
+
+When fitting the data and in the `log_likelihood_function` below, the PSF is used to create the model data. This 
+demonstrates how an `Analysis` class can be extended to include additional steps in the model fitting process.
+"""
+psf = np.load(file=path.join(dataset_path, "psf.npy"))
+plot_array(array=psf, title="Point Spread Function of Galaxy ?")
+
+"""
+__Grid__
+
+To perform certain calculations, we need a grid of (x,y) coordinates that overlap the observed galaxy data.
+
+We create a 2D grid of coordinates where the origin is (0.0, 0.0) and each pixel is 0.05 in size, matching the 
+resolution of our image data.
+
+This grid includes only (y,x) coordinates within the circular mask applied to the data, as we only need to perform 
+calculations within this masked region.
+"""
+grid = np.load(file=path.join(dataset_path, "grid.npy"))
+
+plot_grid(
+    grid=grid,
+    title="Cartesian grid of (x,y) coordinates aligned with dataset",
+)
+
+"""
+__Light Profiles__
+
+Our galaxy model must describe the light of each galaxy, which we refer to as a "light profile." Below, we define two 
+light profile classes named `LightBulgey` and `LightDisky`.
+
+These Python classes serve as the model components representing each galaxy's light, similar to the `Gaussian` class 
+in previous tutorials. The `__init__` constructor's input parameters (e.g., `centre`, `axis_ratio`, `angle`) are the 
+model parameters that the non-linear search will fit.
+
+These classes also include functions that create an image from the light profile based on an input grid of (x,y) 2D 
+coordinates, which we will use to generate an image of a light profile. There is a lot of maths going on in these
+functions which understanding in detail isn't necessary for you to use **PyAutoFit**. So, if the code below looks
+complicated, don't worry about it!
+
+The Python code below uses Python's class inheritance functionality, which novice users may not be familiar with.
+Inheritance is the part of the code which reads `class LightProfile(GeometryProfile):`. This means that 
+the `LightProfile` class inherits the functions and attributes of the `GeometryProfile` class. 
+
+This is a common object-oriented programming feature in Python, but is not something you need to be familiar with to
+use **PyAutoFit**. If the code below is confusing, don't worry about it for now, it won't impact your ability to use
+**PyAutoFit**.
+"""
+
+
+class GeometryProfile:
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+    ):
+        """
+        Abstract base class for the geometry of a profile representing the light or mass of a galaxy.
+
+        Using the centre, axis-ratio and position angle of the profile this class describes how to convert a
+        (y,x) grid of Cartesian coordinates to the elliptical geometry of the profile.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        """
+
+        self.centre = centre
+        self.axis_ratio = axis_ratio
+        self.angle = angle
+
+    def transformed_to_reference_frame_grid_from(self, grid: np.ndarray):
+        """
+        Transform a grid of (y,x) coordinates to the geometric reference frame of the profile via a translation using
+        its `centre` and a rotation using its `angle`.
+
+        This performs the following steps:
+
+        1) Translate the input (y,x) coordinates from the (0.0, 0.0) origin to the centre of the profile by subtracting
+           the profile centre.
+
+        2) Compute the radial distance of every translated coordinate from the centre.
+
+        3) Rotate the coordinates from the above step counter-clockwise from the positive x-axis by the
+           profile's `angle`.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinate grid in its original reference frame.
+        """
+
+        shifted_grid = grid - self.centre
+        effective_radius = (
+            (shifted_grid[:, :, 0] ** 2.0 + shifted_grid[:, :, 1] ** 2.0)
+        ) ** 0.5
+
+        theta_coordinate_to_profile = np.arctan2(
+            shifted_grid[:, :, 1], shifted_grid[:, :, 0]
+        ) - np.radians(self.angle)
+
+        transformed_grid = np.zeros(grid.shape)
+
+        transformed_grid[:, :, 0] = effective_radius * np.cos(
+            theta_coordinate_to_profile
+        )
+        transformed_grid[:, :, 1] = effective_radius * np.sin(
+            theta_coordinate_to_profile
+        )
+
+        return transformed_grid
+
+    def rotated_grid_from_reference_frame_from(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Rotate a grid of (y,x) coordinates which have been transformed to the elliptical reference frame of a profile
+        back to the original unrotated coordinate frame.
+
+        This performs the following steps:
+
+        1) Rotate the coordinates from the elliptical reference frame back to the original unrotated coordinate frame
+           using the profile's `angle`.
+
+        2) Translate the coordinates from the centre of the profile back to the (0.0, 0.0) origin by adding the
+        profile centre.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the reference frame of an elliptical profile.
+        """
+        cos_angle = np.cos(np.radians(self.angle))
+        sin_angle = np.sin(np.radians(self.angle))
+
+        transformed_grid = np.zeros(grid.shape)
+
+        transformed_grid[:, :, 0] = np.add(
+            np.multiply(grid[:, :, 0], cos_angle),
+            -np.multiply(grid[:, :, 1], sin_angle),
+        )
+        transformed_grid[:, :, 1] = np.add(
+            np.multiply(grid[:, :, 0], sin_angle), np.multiply(grid[:, :, 1], cos_angle)
+        )
+
+        return transformed_grid
+
+    def elliptical_radii_grid_from(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Convert a grid of (y,x) coordinates to a grid of elliptical radii.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates in the reference frame of the elliptical profile.
+        """
+
+        return (
+            (grid[:, :, 1] ** 2.0) + (grid[:, :, 0] / self.axis_ratio) ** 2.0
+        ) ** 0.5
+
+
+class LightProfile(GeometryProfile):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        Abstract base class for a light profile, which describes the emission of a galaxy as a
+        function of radius.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(centre=centre, axis_ratio=axis_ratio, angle=angle)
+
+        self.intensity = intensity
+        self.effective_radius = effective_radius
+
+
+class LightBulgey(LightProfile):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        The De Vaucouleurs light profile often used in Astronomy to represent the bulge of galaxies.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(
+            centre=centre,
+            axis_ratio=axis_ratio,
+            angle=angle,
+            intensity=intensity,
+            effective_radius=effective_radius,
+        )
+
+    def image_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the image of the De Vaucouleurs light profile on a grid of Cartesian (y,x) coordinates, which are
+        first translated to the profile's reference frame.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates where the image is computed.
+        """
+        grid_transformed = self.transformed_to_reference_frame_grid_from(grid=grid)
+        grid_elliptical_radii = self.elliptical_radii_grid_from(grid=grid_transformed)
+
+        return self.intensity * np.exp(
+            -7.66924
+            * ((grid_elliptical_radii / self.effective_radius) ** (1.0 / 7.66924) - 1.0)
+        )
+
+
+class LightDisky(LightProfile):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        axis_ratio: float = 1.0,
+        angle: float = 0.0,
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+    ):
+        """
+        The Exponential light profile often used in Astronomy to represent the disk of galaxies.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) coordinates of the profile centre.
+        axis_ratio
+            The axis-ratio of the ellipse (minor axis / major axis).
+        angle
+            The rotation angle in degrees counter-clockwise from the positive x-axis.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile.
+        """
+
+        super().__init__(
+            centre=centre,
+            axis_ratio=axis_ratio,
+            angle=angle,
+            intensity=intensity,
+            effective_radius=effective_radius,
+        )
+
+    def image_from_grid(self, grid: np.ndarray) -> np.ndarray:
+        """
+        Returns the image of the light profile on a grid of Cartesian (y,x) coordinates.
+
+        Parameters
+        ----------
+        grid
+            The (y, x) coordinates where the image is computed.
+        """
+        grid_transformed = self.transformed_to_reference_frame_grid_from(grid=grid)
+        grid_elliptical_radii = self.elliptical_radii_grid_from(grid=grid_transformed)
+
+        return self.intensity * np.exp(
+            -1.67838
+            * ((grid_elliptical_radii / self.effective_radius) ** (1.0 / 1.67838) - 1.0)
+        )
+
+
+"""
+Here is an example of an image of a `LightDisky` profile, using the `image_from_grid` method and the
+grid aligned with the image of the galaxy.
+"""
+light_profile = LightDisky(
+    centre=(0.01, 0.01), axis_ratio=0.7, angle=45.0, intensity=1.0, effective_radius=2.0
+)
+light_image = light_profile.image_from_grid(grid=grid)
+
+plot_array(array=light_image, title="Image of an Exponential light profile.")
+
+"""
+__Model Data__
+
+To produce the `model_data`, we now convolve the overall image with the Point Spread Function (PSF) of our observations. 
+This blurs the image to simulate the effects of the telescope optics and pixelization used to capture the image.
+"""
+model_data = signal.convolve2d(light_image, psf, mode="same")
+
+plot_array(array=model_data, title="Model Data of the Light Profile.")
+
+"""
+By subtracting the model image from the data, we create a 2D residual map, similar to the residual maps we made in 
+the 1D Gaussian examples but now for 2D imaging data. It's evident that the random model we used here does not fit 
+the galaxy well.
+"""
+residual_map = data - model_data
+
+plot_array(array=residual_map, title="Residual Map of fit")
+
+"""
+Just like in the 1D `Gaussian` fitting examples, we can use the noise map to compute the normalized residuals and 
+chi-squared map for the model.
+"""
+normalized_residual_map = residual_map / noise_map
+chi_squared_map = (normalized_residual_map) ** 2.0
+
+plot_array(
+    array=normalized_residual_map,
+    title="Normalized Residual Map of fit",
+)
+plot_array(array=chi_squared_map, title="Chi Squared Map of fit")
+
+"""
+Finally, we compute the `log_likelihood` of this model, which will be used next to fit the model to the data using 
+a non-linear search.
+"""
+chi_squared = np.sum(chi_squared_map)
+noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))
+
+log_likelihood = -0.5 * (chi_squared + noise_normalization)
+
+print(log_likelihood)
+
+"""
+__Model__
+
+We now define a model where one of the light profiles defined earlier is encapsulated within a `Model` object. 
+This allows us to treat the light profile as a model component that can be used in conjunction with other components. 
+
+Although we are currently using only one light profile in this example, we structure it within a `Collection` for 
+potential extension with multiple light profiles.
+"""
+model = af.Collection(light_0=af.Model(LightDisky))
+
+"""
+The `model` operates the same as the model components we've utilized previously. 
+
+Fitting 2D data is more time-consuming than 1D Gaussian data. I therefore employ the prior tuning method described in 
+tutorial 4 to speed up the fit, by adjusting the priors to approximate values close to the correct solution. 
+This adjustment is based on thorough fits that were performed using a slow and detailed search in order to locate
+the global maxima of the likelihood.
+
+Additionally, I fix the center of the light profile to (0.0, 0.0), where it visually appears to be located. 
+This reduces the number of free parameters and simplifies the complexity of the non-linear search.
+"""
+model.light_0.centre_0 = 0.0
+model.light_0.centre_1 = 0.0
+model.light_0.angle = af.UniformPrior(lower_limit=100.0, upper_limit=120.0)
+model.light_0.axis_ratio = af.UniformPrior(lower_limit=0.6, upper_limit=0.8)
+model.light_0.effective_radius = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)
+
+"""
+The model info contains information on all of the model components and priors, including the updates above.
+"""
+print(model.info)
+
+"""
+__Analysis__
+
+We now define the `Analysis` class for this astronomy example, which will fit the `model` to the image data.
+
+Checkout all the docstrings in this class for details on how the fit is performed to 2D imaging data, including
+the role of the mask and PSF.
+"""
+
+
+class Analysis(af.Analysis):
+    def __init__(
+        self,
+        data: np.ndarray,
+        noise_map: np.ndarray,
+        psf: np.ndarray,
+        grid: np.ndarray,
+        mask: np.ndarray,
+    ):
+        """
+        The analysis class for the **PyAutoFit** example Astronomy project on galaxy fitting.
+
+        This class contains imaging data of a galaxy and it fits it with a model which represents the light profile of
+        the galaxy.
+
+        The imaging data contains:
+
+        1) An image of the galaxy.
+
+        2) The noise in every pixel of that image.
+
+        3) The Point Spread Function (PSF) describing how the optics of the telescope blur the image.
+
+        4) A (y,x) grid of coordinates describing the locations of the image pixels in a unit system, which are used
+           to perform calculations.
+
+        5) A mask that removes certain image pixels from the analysis.
+
+        This project is a scaled down version of the Astronomy project **PyAutoGalaxy**, which is one of the
+        original projects from which PyAutoFit is an offshoot!
+
+        https://github.com/Jammy2211/PyAutoGalaxy
+
+        Parameters
+        ----------
+        data
+            The image containing the observation of the galaxy that is fitted.
+        noise_map
+            The RMS noise values of the image data, which is folded into the log likelihood calculation.
+        psf
+            The Point Spread Function of the telescope, which describes how the telescope blurs the image.
+        grid
+            The (y, x) coordinates of the image from which the calculation is performed and model image is
+            computed using.
+        mask
+            The 2D mask that is applied to the image data.
+        """
+
+        super().__init__()
+
+        self.data = data
+        self.noise_map = noise_map
+        self.psf = psf
+        self.grid = grid
+        self.mask = mask
+
+    def log_likelihood_function(self, instance) -> float:
+        """
+        The `log_likelihood_function` of the galaxy example, which performs the following step:
+
+        1) Using the model passed into the function (whose parameters are set via the non-linear search) create
+        a model image of this galaxy using the light profiles in the `instance`.
+
+        2) Convolve this model image with the Point Spread Function of the telescope, ensuring the telescope optics
+        are included in the model image.
+
+        3) Subtract this model image from the data and compute its residuals, chi-squared and likelihood via the
+        noise map.
+
+        4) Apply the mask to all these quantities, ensuring the edges of the galaxy are omitted from the fit.
+
+        Parameters
+        ----------
+        instance
+            An instance of the model set via the non-linear search.
+
+        Returns
+        -------
+        float
+            The log likelihood value of this particular model.
+        """
+
+        """
+        The 'instance' that comes into this method contains the light profiles we setup in the `Model` and `Collection`,
+        which can be seen by uncommenting the code below.
+        """
+
+        # print("Model Instance:")
+        # print("Light Profile = ", instance.light)
+        # print("Light Profile Centre= ", instance.light.centre)
+
+        """
+        Generate the model data from the instance using the functions of the light profile class above.
+        
+        See the docstring of the `model_data_from_instance` function for a description of how this works.
+        """
+
+        model_data = self.model_data_from_instance(instance=instance)
+
+        """
+        In this context of fitting 2D astronomical imaging data, the calculation of residual-map, 
+        normalized residual-map, chi-squared-map, and likelihood differs from previous examples such as fitting 1D 
+        Gaussians. 
+        
+        The main difference lies in the incorporation of a mask, which excludes regions where the galaxy is not 
+        visible, typically the edges.
+
+        To ensure that the mask influences these calculations correctly, we use the `where` parameter of numpy. 
+        This parameter allows us to compute the chi-squared and likelihood only where the mask indicates valid 
+        data (where the mask is 0, indicating unmasked regions). This approach ensures that the analysis focuses 
+        exclusively on the relevant parts of the galaxy image where meaningful comparisons between the model and 
+        data can be made.
+        """
+        residual_map = np.subtract(
+            self.data,
+            model_data,
+            out=np.zeros_like(data),
+            where=np.asarray(self.mask) == 0,
+        )
+
+        normalized_residual_map = np.divide(
+            residual_map,
+            self.noise_map,
+            out=np.zeros_like(residual_map),
+            where=np.asarray(self.mask) == 0,
+        )
+
+        chi_squared_map = np.square(
+            np.divide(
+                residual_map,
+                self.noise_map,
+                out=np.zeros_like(residual_map),
+                where=np.asarray(mask) == 0,
+            )
+        )
+
+        chi_squared = float(np.sum(chi_squared_map[np.asarray(self.mask) == 0]))
+        noise_normalization = float(
+            np.sum(np.log(2 * np.pi * noise_map[np.asarray(mask) == 0] ** 2.0))
+        )
+
+        log_likelihood = -0.5 * (chi_squared + noise_normalization)
+
+        return log_likelihood
+
+    def model_data_from_instance(self, instance):
+        """
+        Create the image of a galaxy, including blurring due to the Point Spread Function of the telescope.
+
+        For the purpose of illustrating **PyAutoFit** you do not need to understand what this function is doing, the
+        main thing to note is that it allows us to create a model image of a galaxy.
+
+        Nevertheless, if you are curious, it inputs the (y,x) Cartesian grids into the light profiles and evaluates
+        the image of the galaxy. Multiple light profiles may be summed to create the overall model image that is fitted.
+
+        The mask is used to zero all values of the model image that are not included in the mask.
+
+        Parameters
+        ----------
+        instance
+            An instance of the model set via the non-linear search.
+
+        Returns
+        -------
+        image
+            The image of this galaxy light profile model.
+        """
+
+        overall_image = np.zeros(self.data.shape)
+
+        for light_profile in instance:
+            overall_image += light_profile.image_from_grid(grid=self.grid)
+
+        model_data = signal.convolve2d(overall_image, self.psf, mode="same")
+
+        model_data[self.mask == 1] = 0.0
+
+        return model_data
+
+
+"""
+__Model Fit__
+
+We have successfully composed a model and analysis class for the astronomy example. 
+
+We can now fit it to the data using a search. We use the same method as previous examples, nested sampling 
+algorithm Dynesty.
+"""
+search = af.DynestyStatic(
+    nlive=500,
+    sample="rwalk",  # This makes dynesty run faster, don't worry about what it means for now!
+)
+
+analysis = Analysis(data=data, noise_map=noise_map, psf=psf, grid=grid, mask=mask)
+
+print(
+    """
+    The non-linear search has begun running.
+    This Jupyter notebook cell with progress once the search has completed - this could take a few minutes!
+    """
+)
+
+result = search.fit(model=model, analysis=analysis)
+
+print("The search has finished run - you may now continue the notebook.")
+
+
+"""
+__Result__
+
+The `result` object provides a concise summary of the outcomes from the non-linear search. While many of these values 
+may not be directly interpretable for non-astronomers, some are meaningful even without specialized knowledge.
+
+For instance, the `angle` parameter indicates the rotation of the galaxy's light profile counterclockwise from the 
+x-axis. In our fit, this value should approximate 110.0 degrees, aligning with the observed orientation of the 
+galaxy in the image.
+
+Similarly, the `axis-ratio` represents the ratio of the minor axis to the major axis of the galaxy's light profile. 
+A value near 0.7 is expected, reflecting the elliptical shape typical of galaxies.
+
+This illustrates why we perform model-fitting, we can take complex data and infer simple, interpretable properties
+from it which provide insight into the physical processes generating the data. This is the core goal of the scientific
+method and the use of models to explain observations.
+"""
+print(result.info)
+
+"""
+The best way to assess the quality of the fit is to use the visualization techniques introduced in previous tutorials. 
+
+For 2D data, we visualize the model image, residual-map, normalized residual-map, and chi-squared-map as 2D arrays. 
+These quantities still have the same meanings as their 1D counterparts. For example, normalized residual-map values 
+still represent the sigma values of the residuals, with values greater than 3.0 indicating a poor fit to the data.
+"""
+model_data = analysis.model_data_from_instance(
+    instance=result.max_log_likelihood_instance
+)
+
+residual_map = np.subtract(
+    data, model_data, out=np.zeros_like(data), where=np.asarray(mask) == 0
+)
+
+normalized_residual_map = np.divide(
+    residual_map,
+    noise_map,
+    out=np.zeros_like(residual_map),
+    where=np.asarray(mask) == 0,
+)
+
+chi_squared_map = np.square(
+    np.divide(
+        residual_map,
+        noise_map,
+        out=np.zeros_like(residual_map),
+        where=np.asarray(mask) == 0,
+    )
+)
+
+plot_array(
+    array=model_data,
+    title="Model Data of the Light Profile.",
+)
+plot_array(
+    array=residual_map,
+    title="Residual Map of fit",
+)
+plot_array(
+    array=normalized_residual_map,
+    title="Normalized Residual Map of fit",
+)
+plot_array(
+    array=chi_squared_map,
+    title="Chi Squared Map of fit",
+)
+
+"""
+__Bulgey__
+
+The fit above utilized the disky light profile, which is typically suitable for disk-like late-type galaxies. 
+
+Now, we will repeat the fit using the bulgey light profile, which is more suitable for bulge-like early-type galaxies.
+
+The fit with the higher log likelihood will provide insight into whether the galaxy is more likely to be an early-type 
+or a late-type galaxy.
+"""
+result_disk = result
+
+model = af.Collection(light_0=af.Model(LightBulgey))
+
+# model.light_0.centre_0 = 0.0
+# model.light_0.centre_1 = 0.0
+# model.light_0.axis_ratio = af.UniformPrior(lower_limit=0.7, upper_limit=1.0)
+
+print(
+    """
+    The non-linear search has begun running.
+    This Jupyter notebook cell with progress once the search has completed - this could take a few minutes!
+    """
+)
+
+result = search.fit(model=model, analysis=analysis)
+
+print("The search has finished run - you may now continue the notebook.")
+
+"""
+Print the result info of the bulgey fit.
+"""
+print(result.info)
+
+"""
+We perform the same visualization of the model image, residual-map, normalized residual-map, and chi-squared-map as
+before.
+
+The model image of the bulgey profile provides a better fit to the data than the disky profile, as the residuals are
+lower and the chi-squared-map values are closer to 0.0.
+
+This suggests that the galaxy is more likely to be an early-type galaxy with a bulge-like light profile.
+"""
+
+model_data = analysis.model_data_from_instance(
+    instance=result.max_log_likelihood_instance
+)
+
+residual_map = np.subtract(
+    data, model_data, out=np.zeros_like(data), where=np.asarray(mask) == 0
+)
+
+normalized_residual_map = np.divide(
+    residual_map,
+    noise_map,
+    out=np.zeros_like(residual_map),
+    where=np.asarray(mask) == 0,
+)
+
+chi_squared_map = np.square(
+    np.divide(
+        residual_map,
+        noise_map,
+        out=np.zeros_like(residual_map),
+        where=np.asarray(mask) == 0,
+    )
+)
+
+plot_array(
+    array=model_data,
+    title="Model Data of the Light Profile.",
+)
+plot_array(
+    array=residual_map,
+    title="Residual Map of fit",
+)
+plot_array(
+    array=normalized_residual_map,
+    title="Normalized Residual Map of fit",
+)
+plot_array(
+    array=chi_squared_map,
+    title="Chi Squared Map of fit",
+)
+
+"""
+To make certain of our interpretation, we should compare the log likelihoods of the two fits. 
+
+The fit with the highest log likelihood is the preferred model, which (provided your non-linear search sampled 
+parameter space accurately), is the bulgey profile.
+
+Therefore, the galaxy is likely an early-type galaxy with a bulge-like light profile.
+"""
+print("Disk Model Log Likelihood:")
+print(result_disk.log_likelihood)
+print("Bulge Model Log Likelihood:")
+print(result.log_likelihood)
+
+"""
+__Model Mismatch__
+
+The analysis above allowed us to determine whether the galaxy is more likely to be an early-type or late-type galaxy.
+
+However, after fitting the bulgey profile, you may not of expected it to be the highest log likelihood model. The
+model gave a relatively poor fit, with significant residuals and chi-squared values. It just turns out that
+the disky profile gave an even worse fit!
+
+This reflected the notion of "model mismatch" that we discussed in tutorial 4. One of the challenges of model fitting
+is that you may not have a model that is a brilliant representation of the data, and your search is successfully
+locating the global maxima even though the fit looks visibly poor.
+
+In Astronomy, what a scientist would do next is update their model to try and improve the fit. For example, they 
+may extend the model to contain both the bulgey and disky profile, allowing the model to fit both components of the
+galaxy simultaneously. 
+
+There are a whole range of approaches that Astronomers take to improve the model, which include fitting the
+galaxy with hundreds of 2D Gaussians, fitting it with a pixel grid of light and decomposing it into what
+are called basis functions. These approaches are not relevent to your understanding of **PyAutoFit**, but they
+are all implemented in **PyAutoGalaxy** via the **PyAutoFit** API for model composition. 
+
+This is another great example of how **PyAutoFit** can be used to perform complex model-fitting tasks.
+
+__Extensions__
+
+To conclude, I will illustrate some of the extensions that can be made to this example and the challenges that arise
+when fitting more complex models.
+
+I don't provide code examples of how to tackle these extensions with model-fitting, but I encourage you to think about
+how **PyAutoFit** could be used to address these problems and go ahead and give them a go yourself if you're feeling
+adventurous!
+
+**Multiple Components**: 
+
+The model above fitted a single light profile to the galaxy, which turned out to be a bulge-like component. 
+
+However, galaxies often have multiple components, for example a disk and bulge, but also structures such as a bar or
+spiral arms. Here is a galaxy with a bulge, disk, bar and spiral arms:
+
+![ngc](https://github.com/Jammy2211/autofit_workspace/blob/main/scripts/howtofit/chapter_1_introduction/images/ngc1300.jpg?raw=true)
+ 
+To model this galaxy, we would add many different Python classes with their own parameters and light profiles 
+representing each component of the galaxy. They would then be combined into a single model, whose model image
+would be the summed image of each component, the total number of parameters would be in the twenties or thirties
+and the degeneracies between the parameters would be challenging to sample accurately. 
+
+Here is a snippet of rougyly what our model composition would look like:
+
+model = af.Collection(
+    bulge=LightBulge,
+    disk=LightDisk,
+    bar=LightBar,
+    spiral_arms=LightSpiralArms
+)
+
+This is a great example of how quickly model-fitting can become complex and how keeping track of the speed and
+efficiency of the non-linear search becomes crucial. Furthermore, with so many different light profiles to fit,
+**PyAutoFit**'s model composition API becomes invaluable.
+
+
+**Multiple Galaxies**:
+
+There is no reason our imaging data need contain only one galaxy. The data could include multiple galaxies, each
+of which we want to fit with its own light profile.
+
+Two galaxies close together are called galaxy mergers, and here is a beautiful example of a pair of merging galaxies:
+
+![merger](https://github.com/Jammy2211/autofit_workspace/blob/main/scripts/howtofit/chapter_1_introduction/images/merger.jpg?raw=true)
+
+This would pretty much double the model complexity, as each galaxy would have its own model and light profile.
+All the usual issues then become doubly important, such as ensuring the likelihood function is efficient, that
+the search can sample parameter space accurately and that the results are interpreted correctly.
+
+The model composition would again change, and might look something like:
+
+model = af.Collection(
+    bulge_0=LightBulge,
+    disk_0=LightDisk,
+    bulge_1=LightBulge,
+    disk_1=LightDisk
+)
+
+model.bulge_0.centre_0 = 2.0
+model.bulge_0.centre_1 = 0.0
+model.disk_0.centre_0 = 2.0
+model.disk_0.centre_1 = 0.0
+
+model.bulge_1.centre_0 = -2.0
+model.bulge_1.centre_1 = 0.0
+model.disk_1.centre_0 = -2.0
+model.disk_1.centre_1 = 0.0
+
+In order to keep the model slightly more simple, the centres of the light profiles have been fixed to where
+they peak in the image. 
+
+**PyAutoFit** extensible model composition API actually has much better tools for composing complex models like this
+than the example above. You can find a concise run through of these in the model cookbook, but they will
+also be the focus on a tutorial in the next chapter of the **HowToFit** lectures.
+
+
+**Multiple Wavelengths**:
+
+Galaxies emit light in many wavelengths, for example ultraviolet, optical, infrared, radio and X-ray. Each wavelength
+provides different information about the galaxy, for example the ultraviolet light tells us about star formation,
+the optical light about the stars themselves and the infrared about dust.
+
+The image below shows observations of the same galaxy at different wavelengths:
+
+![multiband](https://github.com/Jammy2211/autofit_workspace/blob/main/scripts/howtofit/chapter_1_introduction/images/overview_3.jpg?raw=true)
+
+In tutorial 6, we learn how to perform fits to multiple datasets simultaneously, and how to change the model
+parameterization to account for the variation across datasets. 
+
+Multi-wavelength modeling of galaxies is a great example of where this is useful, as it allows us to fit the galaxy
+with certain parameters shared across all wavelengths (e.g., its centre) and other parameters varied (e.g., its
+intensity and effective radius). A child project of **PyAutoFit**, called **PyAutoGalaxy**, uses this functionality
+to achieve exactly this.
+
+__Chapter Wrap Up__
+
+We have now completed the first chapter of **HowToFit**, which has taught you the basics of model-fitting.
+
+Its now time you take everything you've learnt and apply it to your own model-fitting problem. Think carefully
+about the key concepts of this chapter, for example how to compose a model, how to create an analysis class and
+how to overcome challenges that arise when fitting complex models.
+
+Once you are familiar with these concepts, and confident that you have some simple model-fitting problems under
+your belt, you should move on to the next chapter. This covers how to build a scientific workflow around your
+model-fitting, so that you can begin to scale up your model-fitting to more complex models, larger datasets and
+more difficult problems. Checkout the `start_here.ipynb` notebook in chapter 2 to get started!
+"""