notebooks/imaging/advanced/log_likelihood_function/parametric.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Log Likelihood Function: Inversion (Parametric)__\n",
        "\n",
        "This script provides a step-by-step guide of the **PyAutoLens** `log_likelihood_function` which is used to fit\n",
        "`Imaging` data with a parametric lens light profile and source light profile (specifically \n",
        "an `Sersic`, elliptical Sersic, lens and source).\n",
        "\n",
        "This script has the following aims:\n",
        "\n",
        " - To provide a resource that authors can include in papers using **PyAutoLens**, so that readers can understand the\n",
        " likelihood function (including references to the previous literature from which it is defined) without having to\n",
        " write large quantities of text and equations.\n",
        "\n",
        "Accompanying this script is the `contributor_guide.py` which provides URL's to every part of the source-code that\n",
        "is illustrated in this guide. This gives contributors a linear run through of what source-code functions, modules and\n",
        "packages are called when the likelihood is evaluated."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%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 matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from os import path\n",
        "\n",
        "import autolens as al\n",
        "import autolens.plot as aplt"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Dataset__\n",
        "\n",
        "In order to perform a likelihood evaluation, we first load a dataset.\n",
        "\n",
        "This example fits a simulated strong lens which is simulated using a 0.1 arcsecond-per-pixel resolution (this is lower\n",
        "resolution than the best quality Hubble Space Telescope imaging and close to that of the Euclid space satellite)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_path = path.join(\"dataset\", \"imaging\", \"simple\")\n",
        "\n",
        "dataset = al.Imaging.from_fits(\n",
        "    data_path=path.join(dataset_path, \"data.fits\"),\n",
        "    psf_path=path.join(dataset_path, \"psf.fits\"),\n",
        "    noise_map_path=path.join(dataset_path, \"noise_map.fits\"),\n",
        "    pixel_scales=0.1,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Throughout this guide, I will use **PyAutoLens**'s in-built visualization tools for plotting. \n",
        "\n",
        "For example, using the `ImagingPlotter` I can plot the imaging dataset we performed a likelihood evaluation on."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_plotter = aplt.ImagingPlotter(dataset=dataset)\n",
        "dataset_plotter.subplot_dataset()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Mask__\n",
        "\n",
        "The likelihood is only evaluated using image pixels contained within a 2D mask, which we choose before performing\n",
        "lens modeling.\n",
        "\n",
        "Below, we define a 2D circular mask with a 3.0\" radius."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "mask = al.Mask2D.circular(\n",
        "    shape_native=dataset.shape_native, pixel_scales=dataset.pixel_scales, radius=3.0\n",
        ")\n",
        "\n",
        "masked_dataset = dataset.apply_mask(mask=mask)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "When we plot the masked imaging, only the circular masked region is shown."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_plotter = aplt.ImagingPlotter(dataset=masked_dataset)\n",
        "dataset_plotter.subplot_dataset()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Over Sampling__\n",
        "\n",
        "Over sampling evaluates a light profile using multiple samples of its intensity per image-pixel.\n",
        "\n",
        "For simplicity, we disable over sampling in this guide by setting `sub_size=1`. \n",
        "\n",
        "a full description of over sampling and how to use it is given in `autogalaxy_workspace/*/guides/over_sampling.py`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "masked_dataset = masked_dataset.apply_over_sampling(\n",
        "    over_sampling=al.OverSamplingUniform(sub_size=1)\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Masked Image Grid__\n",
        "\n",
        "To perform lensing calculations we first must define the 2D image-plane (y,x) coordinates used in the calculation.\n",
        "\n",
        "These are given by `masked_dataset.grid`, which we can plot and see is a uniform grid of (y,x) Cartesian coordinates\n",
        "which have had the 3.0\" circular mask applied."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "grid_plotter = aplt.Grid2DPlotter(grid=masked_dataset.grid)\n",
        "grid_plotter.figure_2d()\n",
        "\n",
        "print(\n",
        "    f\"(y,x) coordinates of first ten unmasked image-pixels {masked_dataset.grid[0:9]}\"\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To perform lensing calculations we convert this 2D (y,x) grid of coordinates to elliptical coordinates:\n",
        "\n",
        " $\\eta = \\sqrt{(x - x_c)^2 + (y - y_c)^2/q^2}$\n",
        "\n",
        "Where:\n",
        "\n",
        " - $y$ and $x$ are the (y,x) arc-second coordinates of each unmasked image-pixel, given by `masked_dataset.grid`.\n",
        " - $y_c$ and $x_c$ are the (y,x) arc-second `centre` of the light or mass profile used to perform lensing calculations.\n",
        " - $q$ is the axis-ratio of the elliptical light or mass profile (`axis_ratio=1.0` for spherical profiles).\n",
        " - The elliptical coordinates is rotated by position angle $\\phi$, defined counter-clockwise from the positive \n",
        " x-axis.\n",
        "\n",
        "**PyAutoLens** does not use $q$ and $\\phi$ to parameterize the lens model but expresses these as `elliptical_components`:\n",
        "\n",
        "$\\epsilon_{1} =\\frac{1-q}{1+q} \\sin 2\\phi, \\,\\,$\n",
        "$\\epsilon_{2} =\\frac{1-q}{1+q} \\cos 2\\phi.$\n",
        "\n",
        "Note that `Ell` is used as shorthand for elliptical and `Sph` for spherical."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "profile = al.EllProfile(centre=(0.1, 0.2), ell_comps=(0.1, 0.2))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we transform `masked_dataset.grid ` to the centre of profile and rotate it using its angle `phi`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "transformed_grid = profile.transformed_to_reference_frame_grid_from(\n",
        "    grid=masked_dataset.grid\n",
        ")\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=transformed_grid)\n",
        "grid_plotter.figure_2d()\n",
        "print(\n",
        "    f\"transformed coordinates of first ten unmasked image-pixels {transformed_grid[0:9]}\"\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using these transformed (y',x') values we compute the elliptical coordinates $\\eta = \\sqrt{(x')^2 + (y')^2/q^2}$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "elliptical_radii = profile.elliptical_radii_grid_from(grid=transformed_grid)\n",
        "\n",
        "print(\n",
        "    f\"elliptical coordinates of first ten unmasked image-pixels {elliptical_radii[0:9]}\"\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Setup: Lens Galaxy Light (Setup)__\n",
        "\n",
        "To perform a likelihood evaluation we now compose our lens model.\n",
        "\n",
        "We first define the light profiles which represents the lens galaxy's light, which will be used to fit the lens \n",
        "light.\n",
        "\n",
        "A light profile is defined by its intensity $I (\\eta_{\\rm l}) $, for example the Sersic profile:\n",
        "\n",
        "$I_{\\rm  Ser} (\\eta_{\\rm l}) = I \\exp \\bigg\\{ -k \\bigg[ \\bigg( \\frac{\\eta}{R} \\bigg)^{\\frac{1}{n}} - 1 \\bigg] \\bigg\\}$\n",
        "\n",
        "Where:\n",
        "\n",
        " - $\\eta$ are the elliptical coordinates (see above) or the masked image-grid.\n",
        " - $I$ is the `intensity`, which controls the overall brightness of the Sersic profile.\n",
        " - $n$ is the ``sersic_index``, which via $k$ controls the steepness of the inner profile.\n",
        " - $R$ is the `effective_radius`, which defines the arc-second radius of a circle containing half the light.\n",
        "\n",
        "In this example, we assume our lens is composed of two light profiles, an elliptical Sersic and Exponential (a Sersic\n",
        "where `sersic_index=4`) which represent the bulge of the lens. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "bulge = al.lp.Sersic(\n",
        "    centre=(0.0, 0.0),\n",
        "    ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=45.0),\n",
        "    intensity=4.0,\n",
        "    effective_radius=0.6,\n",
        "    sersic_index=3.0,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the masked 2D grid defined above, we can calculate and plot images of each light profile component.\n",
        "\n",
        "(The transformation to elliptical coordinates above are built into the `image_2d_from` function and performed \n",
        "implicitly)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "image_2d_bulge = bulge.image_2d_from(grid=masked_dataset.grid)\n",
        "\n",
        "bulge_plotter = aplt.LightProfilePlotter(light_profile=bulge, grid=masked_dataset.grid)\n",
        "bulge_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Setup: Lens Galaxy Mass__\n",
        "\n",
        "We next define the mass profiles which represents the lens galaxy's mass, which will be used to ray-trace the \n",
        "image-plane 2D grid of (y,x) coordinates to the source-plane so that the source model can be evaluated.\n",
        "\n",
        "In this example, we assume our lens is composed of an elliptical isothermal mass distribution and external shear.\n",
        "\n",
        "A mass profile is defined by its convergence $\\kappa (\\eta)$, which is related to\n",
        "the surface density of the mass distribution as\n",
        "\n",
        "$\\kappa(\\eta)=\\frac{\\Sigma(\\eta)}{\\Sigma_\\mathrm{crit}},$\n",
        "\n",
        "where\n",
        "\n",
        "$\\Sigma_\\mathrm{crit}=\\frac{{\\rm c}^2}{4{\\rm \\pi} {\\rm G}}\\frac{D_{\\rm s}}{D_{\\rm l} D_{\\rm ls}},$\n",
        "\n",
        "and\n",
        "\n",
        " - `c` is the speed of light.\n",
        " - $D_{\\rm l}$, $D_{\\rm s}$, and $D_{\\rm ls}$ are respectively the angular diameter distances to the lens, to the \n",
        " source, and from the lens to the source.\n",
        "\n",
        "For readers less familiar with lensing, we can think of $\\kappa(\\eta)$ as a convenient and\n",
        "dimensionless way to describe how light is gravitationally lensed after assuming a cosmology.\n",
        "\n",
        "For the for the isothermal profile:\n",
        "\n",
        "$\\kappa(\\eta) = \\frac{1.0}{1 + q} \\bigg( \\frac{\\theta_{\\rm E}}{\\eta} \\bigg)$\n",
        "\n",
        "Where:\n",
        "\n",
        " - $\\theta_{\\rm E}$ is the `einstein_radius` (which is rescaled compared to other einstein radius\n",
        " definitions)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "mass = al.mp.Isothermal(\n",
        "    centre=(0.0, 0.0),\n",
        "    einstein_radius=1.6,\n",
        "    ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=45.0),\n",
        ")\n",
        "\n",
        "shear = al.mp.ExternalShear(gamma_1=0.05, gamma_2=0.05)\n",
        "\n",
        "mass_plotter = aplt.MassProfilePlotter(mass_profile=mass, grid=masked_dataset.grid)\n",
        "mass_plotter.figures_2d(convergence=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From each mass profile we can compute its deflection angles, which describe how due to gravitational lensing\n",
        "image-pixels are ray-traced to the source plane.\n",
        "\n",
        "The deflection angles are computed by integrating $\\kappa$: \n",
        "\n",
        "$\\vec{{\\alpha}}_{\\rm x,y} (\\vec{x}) = \\frac{1}{\\pi} \\int \\frac{\\vec{x} - \\vec{x'}}{\\left | \\vec{x} - \\vec{x'} \\right |^2} \\kappa(\\vec{x'}) d\\vec{x'} \\, ,$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "deflections_yx_2d = mass.deflections_yx_2d_from(grid=masked_dataset.grid)\n",
        "\n",
        "mass_plotter = aplt.MassProfilePlotter(mass_profile=mass, grid=masked_dataset.grid)\n",
        "mass_plotter.figures_2d(deflections_y=True, deflections_x=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Setup: Lens Galaxy__\n",
        "\n",
        "We now combine the light and mass profiles into a single `Galaxy` object for the lens galaxy.\n",
        "\n",
        "When computing quantities for the light and mass profiles from this object, it computes each individual quantity and \n",
        "adds them together. \n",
        "\n",
        "For example, for the `bulge`, when it computes their 2D images it computes each individually and then adds\n",
        "them together."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_galaxy = al.Galaxy(redshift=0.5, bulge=bulge, mass=mass, shear=shear)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Setup: Source Galaxy Light Profile__\n",
        "\n",
        "The source galaxy is fitted using another analytic light profile, in this example another elliptical Sersic."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "source_galaxy = al.Galaxy(\n",
        "    redshift=1.0,\n",
        "    bulge=al.lp.Sersic(\n",
        "        centre=(0.0, 0.0),\n",
        "        ell_comps=al.convert.ell_comps_from(axis_ratio=0.8, angle=60.0),\n",
        "        intensity=4.0,\n",
        "        effective_radius=0.1,\n",
        "        sersic_index=1.0,\n",
        "    ),\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 1: Lens Light__\n",
        "\n",
        "Compute a 2D image of the lens galaxy's light as the sum of its individual light profiles (the `Sersic` \n",
        "bulge). \n",
        "\n",
        "This computes the `lens_image_2d` of each `LightProfile` and adds them together. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_image_2d = lens_galaxy.image_2d_from(grid=masked_dataset.grid)\n",
        "\n",
        "galaxy_plotter = aplt.GalaxyPlotter(galaxy=lens_galaxy, grid=masked_dataset.grid)\n",
        "galaxy_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To convolve the lens's 2D image with the imaging data's PSF, we need its `blurring_image`. This represents all flux \n",
        "values not within the mask, which are close enough to it that their flux blurs into the mask after PSF convolution.\n",
        "\n",
        "To compute this, a `blurring_mask` and `blurring_grid` are used, corresponding to these pixels near the edge of the \n",
        "actual mask whose light blurs into the image:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_blurring_image_2d = lens_galaxy.image_2d_from(grid=masked_dataset.blurring_grid)\n",
        "\n",
        "galaxy_plotter = aplt.GalaxyPlotter(\n",
        "    galaxy=lens_galaxy, grid=masked_dataset.blurring_grid\n",
        ")\n",
        "galaxy_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 2: Ray Tracing__\n",
        "\n",
        "To perform lensing calculations we ray-trace every 2d (y,x) coordinate $\\theta$ from the image-plane to its (y,x) \n",
        "source-plane coordinate $\\beta$ using the summed deflection angles $\\alpha$ of the mass profiles:\n",
        "\n",
        " $\\beta = \\theta - \\alpha(\\theta)$\n",
        "\n",
        "The likelihood function of a source light profile ray-traces two grids from the image-plane to the source-plane:\n",
        "\n",
        " 1) A 2D grid of (y,x) coordinates aligned with the imaging data's image-pixels.\n",
        " \n",
        " 2) The 2D blurring grid (used for the lens light above) which accounts for pixels at the edge of the mask whose\n",
        " light blurs into the mask.\n",
        " \n",
        "The function below computes the 2D deflection angles of the tracer's lens galaxies and subtracts them from the \n",
        "image-plane 2D (y,x) coordinates $\\theta$ of each grid, thus ray-tracing their coordinates to the source plane to \n",
        "compute their $\\beta$ values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])\n",
        "\n",
        "# A list of every grid (e.g. image-plane, source-plane) however we only need the source plane grid with index -1.\n",
        "traced_grid = tracer.traced_grid_2d_list_from(grid=masked_dataset.grid)[-1]\n",
        "\n",
        "mat_plot = aplt.MatPlot2D(axis=aplt.Axis(extent=[-1.5, 1.5, -1.5, 1.5]))\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=traced_grid, mat_plot_2d=mat_plot)\n",
        "grid_plotter.figure_2d()\n",
        "\n",
        "traced_blurring_grid = tracer.traced_grid_2d_list_from(\n",
        "    grid=masked_dataset.blurring_grid\n",
        ")[-1]\n",
        "\n",
        "mat_plot = aplt.MatPlot2D(axis=aplt.Axis(extent=[-1.5, 1.5, -1.5, 1.5]))\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=traced_blurring_grid, mat_plot_2d=mat_plot)\n",
        "grid_plotter.figure_2d()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 3: Source Image__\n",
        "\n",
        "We pass the traced grid and blurring grid of coordinates to the source galaxy to evaluate its 2D image."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "source_image_2d = source_galaxy.image_2d_from(grid=traced_grid)\n",
        "\n",
        "galaxy_plotter = aplt.GalaxyPlotter(galaxy=lens_galaxy, grid=traced_grid)\n",
        "galaxy_plotter.figures_2d(image=True)\n",
        "\n",
        "source_blurring_image_2d = source_galaxy.image_2d_from(grid=traced_blurring_grid)\n",
        "\n",
        "galaxy_plotter = aplt.GalaxyPlotter(galaxy=source_galaxy, grid=traced_blurring_grid)\n",
        "galaxy_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 4: Lens + Source Light Addition__\n",
        "\n",
        "We add the lens and source galaxy images and blurring together, to create an overall image of the strong lens."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "image = lens_image_2d + source_image_2d\n",
        "\n",
        "array_2d_plotter = aplt.Array2DPlotter(array=image)\n",
        "array_2d_plotter.figure_2d()\n",
        "\n",
        "blurring_image_2d = lens_blurring_image_2d + source_blurring_image_2d\n",
        "\n",
        "array_2d_plotter = aplt.Array2DPlotter(array=blurring_image_2d)\n",
        "array_2d_plotter.figure_2d()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 5: Convolution__\n",
        "\n",
        "Convolve the 2D image of the lens and source above with the PSF in real-space (as opposed to via an FFT) using \n",
        "a `Convolver`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "convolved_image_2d = masked_dataset.convolver.convolve_image(\n",
        "    image=image, blurring_image=blurring_image_2d\n",
        ")\n",
        "\n",
        "array_2d_plotter = aplt.Array2DPlotter(array=convolved_image_2d)\n",
        "array_2d_plotter.figure_2d()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 7: Likelihood Function__\n",
        "\n",
        "We now quantify the goodness-of-fit of our lens and source model.\n",
        "\n",
        "We compute the `log_likelihood` of the fit, which is the value returned by the **PyAutoLens** `log_likelihood_function`.\n",
        "\n",
        "The likelihood function for parametric lens modeling consists of two terms:\n",
        "\n",
        " $-2 \\mathrm{ln} \\, \\epsilon = \\chi^2 + \\sum_{\\rm  j=1}^{J} { \\mathrm{ln}} \\left [2 \\pi (\\sigma_j)^2 \\right]  \\, .$\n",
        "\n",
        "We now explain what each of these terms mean.\n",
        "\n",
        "__Likelihood Step 8: Chi Squared__\n",
        "\n",
        "The first term is a $\\chi^2$ statistic, which is defined above in our merit function as and is computed as follows:\n",
        "\n",
        " - `model_data` = `convolved_image_2d`\n",
        " - `residual_map` = (`data` - `model_data`)\n",
        " - `normalized_residual_map` = (`data` - `model_data`) / `noise_map`\n",
        " - `chi_squared_map` = (`normalized_residuals`) ** 2.0 = ((`data` - `model_data`)**2.0)/(`variances`)\n",
        " - `chi_squared` = sum(`chi_squared_map`)\n",
        "\n",
        "The chi-squared therefore quantifies if our fit to the data is accurate or not. \n",
        "\n",
        "High values of chi-squared indicate that there are many image pixels our model did not produce a good fit to the image \n",
        "for, corresponding to a fit with a lower likelihood."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "model_image = convolved_image_2d\n",
        "\n",
        "residual_map = masked_dataset.data - model_image\n",
        "normalized_residual_map = residual_map / masked_dataset.noise_map\n",
        "chi_squared_map = normalized_residual_map**2.0\n",
        "\n",
        "chi_squared = np.sum(chi_squared_map)\n",
        "\n",
        "print(chi_squared)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `chi_squared_map` indicates which regions of the image we did and did not fit accurately."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "chi_squared_map = al.Array2D(values=chi_squared_map, mask=mask)\n",
        "\n",
        "array_2d_plotter = aplt.Array2DPlotter(array=chi_squared_map)\n",
        "array_2d_plotter.figure_2d()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 9: Noise Normalization Term__\n",
        "\n",
        "Our likelihood function assumes the imaging data consists of independent Gaussian noise in every image pixel.\n",
        "\n",
        "The final term ins the likelihood function is therefore a `noise_normalization` term, which consists of the sum\n",
        "of the log of every noise-map value squared. \n",
        "\n",
        "Given the `noise_map` is fixed, this term does not change during the lens modeling process and has no impact on the \n",
        "model we infer."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "noise_normalization = float(np.sum(np.log(2 * np.pi * masked_dataset.noise_map**2.0)))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Likelihood Step 10: Calculate The Log Likelihood!__\n",
        "\n",
        "We made it!\n",
        "\n",
        "We can now, finally, compute the `log_likelihood` of the lens model, by combining the two terms computed above using\n",
        "the likelihood function defined above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "figure_of_merit = float(-0.5 * (chi_squared + noise_normalization))\n",
        "\n",
        "print(figure_of_merit)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Fit__\n",
        "\n",
        "This 11 step process to perform a likelihood function evaluation is what is performed in the `FitImaging` object, which\n",
        "those of you familiar will have seen before."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fit = al.FitImaging(dataset=masked_dataset, tracer=tracer)\n",
        "fit_figure_of_merit = fit.figure_of_merit\n",
        "print(fit_figure_of_merit)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Lens Modeling__\n",
        "\n",
        "To fit a lens model to data, **PyAutoLens** samples the likelihood function illustrated in this tutorial using a\n",
        "non-linear search algorithm.\n",
        "\n",
        "The default sampler is the nested sampling algorithm `Nautilus` (https://github.com/joshspeagle/Nautilus)\n",
        "but **PyAutoLens** supports multiple MCMC and optimization algorithms. \n",
        "\n",
        "__Wrap Up__\n",
        "\n",
        "We have presented a visual step-by-step guide to the **PyAutoLens** parametric likelihood function, which uses analytic\n",
        "light profiles to fit the lens and source light.\n",
        "\n",
        "There are a number of other inputs features which slightly change the behaviour of this likelihood function, which\n",
        "are described in additional notebooks found in this package. In brief, these describe:\n",
        "\n",
        " - **Sub-gridding**: Oversampling the image grid into a finer grid of sub-pixels, which are all individually \n",
        " ray-traced to the source-plane and used to evaluate the light profile more accurately."
      ]
    },
    {
      "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
}