{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Tutorial 3: Inversions\n", "======================\n", "\n", "In the previous two tutorials, we introduced:\n", "\n", " - `Pixelization`'s: which place a pixel-grid in the source-plane.\n", " - `Mappers`'s: which describe how each source-pixel maps to one or more image pixels.\n", "\n", "However, non of this has actually helped us fit strong lens data or reconstruct the source galaxy. This is the subject\n", "of this tutorial, where the process of reconstructing the source's light on the pixelization is called an `Inversion`." ] }, { "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", "from os import path\n", "import autolens as al\n", "import autolens.plot as aplt" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Initial Setup__\n", "\n", "we'll use the same strong lensing data as the previous tutorial, where:\n", "\n", " - The lens galaxy's light is omitted.\n", " - The lens galaxy's total mass distribution is an `Isothermal` and `ExternalShear`.\n", " - The source galaxy's light is an `Sersic`." ] }, { "cell_type": "code", "metadata": {}, "source": [ "dataset_name = \"simple__no_lens_light\"\n", "dataset_path = path.join(\"dataset\", \"imaging\", dataset_name)\n", "\n", "dataset = al.Imaging.from_fits(\n", " data_path=path.join(dataset_path, \"data.fits\"),\n", " noise_map_path=path.join(dataset_path, \"noise_map.fits\"),\n", " psf_path=path.join(dataset_path, \"psf.fits\"),\n", " pixel_scales=0.1,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets create an annular mask which traces the stongly lensed source's ring of light." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mask = al.Mask2D.circular_annular(\n", " shape_native=dataset.shape_native,\n", " pixel_scales=dataset.pixel_scales,\n", " inner_radius=0.5,\n", " outer_radius=2.8,\n", ")\n", "\n", "visuals = aplt.Visuals2D(mask=mask)\n", "\n", "dataset_plotter = aplt.ImagingPlotter(dataset=dataset, visuals_2d=visuals)\n", "dataset_plotter.figures_2d(data=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now create the masked source-plane grid via the tracer, as we did in the previous tutorial." ] }, { "cell_type": "code", "metadata": {}, "source": [ "dataset = dataset.apply_mask(mask=mask)\n", "\n", "lens_galaxy = al.Galaxy(\n", " redshift=0.5,\n", " 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", " shear=al.mp.ExternalShear(gamma_1=0.05, gamma_2=0.05),\n", ")\n", "\n", "tracer = al.Tracer(galaxies=[lens_galaxy, al.Galaxy(redshift=1.0)])\n", "\n", "source_plane_grid = tracer.traced_grid_2d_list_from(grid=dataset.grid)[1]" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "we again use the rectangular pixelization to create the mapper.\n", "\n", "(Ignore the regularization input below for now, we will cover this in the next tutorial)." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mesh = al.mesh.Rectangular(shape=(25, 25))\n", "\n", "pixelization = al.Pixelization(mesh=mesh)\n", "\n", "mapper_grids = pixelization.mapper_grids_from(\n", " mask=mask, source_plane_data_grid=source_plane_grid\n", ")\n", "mapper = al.Mapper(\n", " mapper_grids=mapper_grids,\n", " over_sampler=dataset.over_sampler_pixelization,\n", " regularization=al.reg.Constant(coefficient=1.0),\n", ")\n", "\n", "include = aplt.Include2D(mask=True, mapper_source_plane_data_grid=True)\n", "\n", "mapper_plotter = aplt.MapperPlotter(mapper=mapper, include_2d=include)\n", "mapper_plotter.subplot_image_and_mapper(image=dataset.data)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Pixelization__\n", "\n", "Finally, we can now use the `Mapper` to reconstruct the source via an `Inversion`. I'll explain how this works in a \n", "second, but lets just go ahead and create the inversion first. " ] }, { "cell_type": "code", "metadata": {}, "source": [ "inversion = al.Inversion(dataset=dataset, linear_obj_list=[mapper])" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inversion has reconstructed the source's light on the rectangular pixel grid, which is called the \n", "`reconstruction`. This source-plane reconstruction can be mapped back to the image-plane to produce the \n", "`mapped_reconstructed_image`." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(inversion.reconstruction)\n", "print(inversion.mapped_reconstructed_image)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both of these can be plotted using an `InversionPlotter`.\n", "\n", "It is possible for an inversion to have multiple `Mapper`'s, therefore for certain figures we specify the index \n", "of the mapper we wish to plot. In this case, because we only have one mapper we specify the index 0." ] }, { "cell_type": "code", "metadata": {}, "source": [ "include = aplt.Include2D(mask=True)\n", "\n", "inversion_plotter = aplt.InversionPlotter(inversion=inversion, include_2d=include)\n", "inversion_plotter.figures_2d(reconstructed_image=True)\n", "inversion_plotter.figures_2d_of_pixelization(pixelization_index=0, reconstruction=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "There we have it, we have successfully reconstructed the source using a rectangular pixel-grid. Whilst this source \n", "was simple (a blob of light in the centre of the source-plane), inversions come into their own when fitting sources \n", "with complex morphologies. \n", "\n", "Lets use an inversion to reconstruct a complex source!" ] }, { "cell_type": "code", "metadata": {}, "source": [ "dataset_name = \"source_complex\"\n", "dataset_path = path.join(\"dataset\", \"imaging\", dataset_name)\n", "\n", "dataset = al.Imaging.from_fits(\n", " data_path=path.join(dataset_path, \"data.fits\"),\n", " noise_map_path=path.join(dataset_path, \"noise_map.fits\"),\n", " psf_path=path.join(dataset_path, \"psf.fits\"),\n", " pixel_scales=0.05,\n", ")\n", "\n", "dataset_plotter = aplt.ImagingPlotter(dataset=dataset)\n", "dataset_plotter.figures_2d(data=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "This code is doing all the same as above (setup the mask, galaxy, tracers, mapper, inversion, etc.)." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mask = al.Mask2D.circular_annular(\n", " shape_native=dataset.shape_native,\n", " pixel_scales=dataset.pixel_scales,\n", " inner_radius=0.1,\n", " outer_radius=3.2,\n", ")\n", "\n", "visuals = aplt.Visuals2D(mask=mask)\n", "\n", "dataset_plotter = aplt.ImagingPlotter(dataset=dataset, visuals_2d=visuals)\n", "dataset_plotter.figures_2d(data=True)\n", "\n", "dataset = dataset.apply_mask(mask=mask)\n", "\n", "lens_galaxy = al.Galaxy(\n", " redshift=0.5,\n", " mass=al.mp.Isothermal(\n", " centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)\n", " ),\n", ")\n", "\n", "tracer = al.Tracer(galaxies=[lens_galaxy, al.Galaxy(redshift=1.0)])\n", "\n", "source_plane_grid = tracer.traced_grid_2d_list_from(grid=dataset.grid)[1]\n", "\n", "mapper_grids = mesh.mapper_grids_from(\n", " mask=mask, source_plane_data_grid=source_plane_grid\n", ")\n", "mapper = al.Mapper(\n", " mapper_grids=mapper_grids,\n", " over_sampler=dataset.over_sampler_pixelization,\n", " regularization=al.reg.Constant(coefficient=1.0),\n", ")\n", "\n", "\n", "inversion = al.Inversion(dataset=dataset, linear_obj_list=[mapper])" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets plot the complex source reconstruction." ] }, { "cell_type": "code", "metadata": {}, "source": [ "inversion_plotter = aplt.InversionPlotter(inversion=inversion, include_2d=include)\n", "inversion_plotter.figures_2d(reconstructed_image=True)\n", "inversion_plotter.figures_2d_of_pixelization(pixelization_index=0, reconstruction=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pretty great, huh? If you ran the complex source pipeline in chapter 3, you'll remember that getting a model image \n", "that looked this good simply *was not possible*. With an inversion, we can do this with ease and without having to \n", "perform model-fitting with 20+ parameters for the source's light!\n", "\n", "We will now briefly discuss how an inversion actually works, however the explanation I give in this tutorial will be \n", "overly-simplified. To be good at lens modeling you do not need to understand the details of how an inversion works, you \n", "simply need to be able to use an inversion to model a strong lens. \n", "\n", "To begin, lets consider some random mappings between our mapper`s source-pixels and the image." ] }, { "cell_type": "code", "metadata": {}, "source": [ "visuals = aplt.Visuals2D(pix_indexes=[[445], [285], [313], [132], [11]])\n", "\n", "mapper_plotter = aplt.MapperPlotter(\n", " mapper=mapper, visuals_2d=visuals, include_2d=include\n", ")\n", "mapper_plotter.subplot_image_and_mapper(image=dataset.data)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "These mappings are known before the inversion reconstructs the source galaxy, which means before this inversion is\n", "performed we know two key pieces of information:\n", "\n", " 1) The mappings between every source-pixel and sets of image-pixels.\n", " 2) The flux values in every observed image-pixel, which are the values we want to fit successfully.\n", "\n", "It turns out that with these two pieces of information we can linearly solve for the set of source-pixel fluxes that \n", "best-fit (e.g. maximize the log likelihood) our observed image. Essentially, we set up the mappings between source and \n", "image pixels as a large matrix and solve for the source-pixel fluxes in an analogous fashion to how you would solve a \n", "set of simultaneous linear equations. This process is called a `linear inversion`.\n", "\n", "There are three more things about a linear inversion that are worth knowing:\n", "\n", " 1) When performing fits using light profiles, we discussed how a `model_image` was generated by convolving the light\n", " profile images with the data's PSF. A similar blurring operation is incorporated into the inversion, such that it \n", " reconstructs a source (and therefore image) which fully accounts for the telescope optics and effect of the PSF.\n", "\n", " 2) You may be familiar with image sub-gridding, which splits each image-pixel into a sub-pixel (if you are not \n", " familiar then feel free to checkout the optional **HowToLens** tutorial on sub-gridding. If a sub-grid is used, it is \n", " the mapping between every sub-pixel and source-pixel that is computed and used to perform the inversion. This prevents \n", " aliasing effects degrading the image reconstruction. By default **PyAutoLens** uses sub-gridding of degree 4x4.\n", "\n", " 3) The inversion`s solution is regularized. But wait, that`s what we'll cover in the next tutorial!\n", "\n", "Finally, let me show you how easy it is to fit an image with an `Inversion` using a `FitImaging` object. Instead of \n", "giving the source galaxy a light profile, we simply pass it a `Pixelization` and regularization, and pass it to a \n", "tracer." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixelization = al.Pixelization(\n", " mesh=al.mesh.Rectangular(shape=(40, 40)),\n", " regularization=al.reg.Constant(coefficient=1.0),\n", ")\n", "\n", "source_galaxy = al.Galaxy(redshift=1.0, pixelization=pixelization)\n", "\n", "tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, like before, we pass the imaging and tracer `FitImaging` object. \n", "\n", "We see some pretty good looking residuals, we must be fitting the lensed source accurately! In fact, we can use the\n", "`subplot_of_planes` method to specifically visualize the inversion and plot the source reconstruction." ] }, { "cell_type": "code", "metadata": {}, "source": [ "fit = al.FitImaging(dataset=dataset, tracer=tracer)\n", "\n", "include = aplt.Include2D(mask=True)\n", "\n", "fit_plotter = aplt.FitImagingPlotter(fit=fit, include_2d=include)\n", "fit_plotter.subplot_fit()\n", "fit_plotter.subplot_of_planes(plane_index=1)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Positive Only Solver__\n", "\n", "All pixelized source reconstructions use a positive-only solver, meaning that every source-pixel is only allowed\n", "to reconstruct positive flux values. This ensures that the source reconstruction is physical and that we don't\n", "reconstruct negative flux values that don't exist in the real source galaxy (a common systematic solution in lens\n", "analysis).\n", "\n", "It may be surprising to hear that this is a feature worth pointing out, but it turns out setting up the linear algebra\n", "to enforce positive reconstructions is difficult to make efficient. A lot of development time went into making this\n", "possible, where a bespoke fast non-negative linear solver was developed to achieve this.\n", "\n", "Other methods in the literature often do not use a positive only solver, and therefore suffer from these \n", "unphysical solutions, which can degrade the results of lens model in general.\n", "\n", "__Wrap Up__\n", "\n", "And, we're done, here are a few questions to get you thinking about inversions:\n", "\n", " 1) The inversion provides the maximum log likelihood solution to the observed image. Is there a problem with seeking \n", " the highest likelihood solution? Is there a risk that we're going to fit other things in the image than just the \n", " lensed source galaxy? What happens if you reduce the `coefficient` of the regularization object above to zero?\n", "\n", " 2) The exterior pixels in the rectangular pixel-grid have no image-pixels in them. However, they are still given a \n", " reconstructed flux. Given these pixels do not map to the data, where is this value coming from?\n", " \n", "__Detailed Explanation__\n", "\n", "If you are interested in a more detailed description of how inversions work, then checkout the file\n", "`autolens_workspace/*/imaging/log_likelihood_function/inversion.ipynb` which gives a visual step-by-step\n", "guide of the process alongside equations and references to literature on the subject." ] }, { "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 }