{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Overview: Pixelizations\n", "-----------------------\n", "\n", "Many strongly lensed source galaxies are complex, and they have asymmetric and irregular morphologies. These\n", "morphologies cannot be well approximated by parametric light profiles like a Sersic, or multiple Sersics. Even\n", "techniques like a multi-Gaussian expansion or shapelets cannot capture the most complex of source morphologies.\n", "\n", "A pixelization reconstructs the source's light using an adaptive pixel-grid, where the solution is regularized using a\n", "prior that forces the solution to have a degree of smoothness.\n", "\n", "This means they can reconstruct more complex source morphologies and are better suited to performing detailed analyses\n", "of a lens galaxy's mass." ] }, { "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": [ "__Dataset__\n", "\n", "We load the imaging data that we'll reconstruct the lensed source galaxy's light of using a pixelization.\n", "\n", "Note how complex the lensed source galaxy looks, with multiple clumps of light. This would be very difficult to \n", "model using light profiles!" ] }, { "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", " psf_path=path.join(dataset_path, \"psf.fits\"),\n", " noise_map_path=path.join(dataset_path, \"noise_map.fits\"),\n", " pixel_scales=0.05,\n", ")\n", "\n", "dataset_plotter = aplt.ImagingPlotter(dataset=dataset)\n", "dataset_plotter.subplot_dataset()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Mask__\n", "\n", "We are going to fit this data via a pixelization, which requires us to define a 2D mask within which the pixelization\n", "reconstructs the data." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mask = al.Mask2D.circular(\n", " shape_native=dataset.shape_native, pixel_scales=dataset.pixel_scales, radius=3.6\n", ")\n", "\n", "dataset = dataset.apply_mask(mask=mask)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Mesh + Regularization__\n", "\n", "To reconstruct the source on a pixel-grid, called a mesh, we simply pass it the `Mesh` class we want to reconstruct its \n", "light. We also pass it a `Regularization` scheme, describing a prior on how smooth the source reconstruction should be. \n", "\n", "Lets use a `Rectangular` mesh with resolution 40 x 40 pixels and `Constant` regularizaton scheme with a \n", "regularization coefficient of 1.0. The higher this coefficient, the more our source reconstruction is smoothed.\n", "\n", "The isothermal mass model defined below is true model used to simulate the data." ] }, { "cell_type": "code", "metadata": {}, "source": [ "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", "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)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Fit__\n", "\n", "Now that our source-galaxy has a `Pixelization`, we are able to fit the data using the same tools described in \n", "a previous overview example. \n", "\n", "We simply pass the source galaxy to a `Tracer` and using this `Tracer` to create a `FitImaging` object." ] }, { "cell_type": "code", "metadata": {}, "source": [ "tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])\n", "\n", "fit = al.FitImaging(dataset=dataset, tracer=tracer)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Pixelization__\n", "\n", "The fit has been performed using a pixelization for the source galaxy, with the following worth noting:\n", "\n", " - The central-right and bottom-right panel shows a pixelized grid of the subplot show the source has been \n", " reconstructed on an uniform rectangular grid of pixels.\n", "\n", " - The source reconstruction is irregular and has multiple clumps of light, these features would be difficult to \n", " represent using analytic light profiles!\n", "\n", " - The source reconstruction has been mapped back to the image-plane, to produce the reconstructed model image, which \n", " is how a `log_likelihood` is computed.\n", " \n", " - This reconstructed model image produces significant residuals, because a rectangular mesh is not an optimal way to\n", " reconstruct the source galaxy." ] }, { "cell_type": "code", "metadata": {}, "source": [ "fit_plotter = aplt.FitImagingPlotter(fit=fit)\n", "fit_plotter.subplot_fit()" ], "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", "__Alternative Pixelizations__\n", "\n", "**PyAutoLens** supports many different meshes. Below, we use a `Delaunay` mesh, which uses a Delaunay tessellation to\n", "reconstruct the source.\n", "\n", "This requires that the centrals of the Delaunay cells, which act as source pixels, are computed first. This uses an\n", "\"image-mesh\", which computes them by overlaying a uniform grid of Cartesian coordinates over the image-plane image\n", "and ray-traces them to the source plane.\n", "\n", "The source pixel-grid is therefore adapted to the mass-model magnification pattern, placing more source-pixel in the\n", "highly magnified regions of the source-plane.\n", "\n", "This leads to a noticeable improvement in the fit, where the residuals are reduced and the source-reconstruction\n", "is noticeably smoother." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixelization = al.Pixelization(\n", " image_mesh=al.image_mesh.Overlay(shape=(40, 40)),\n", " mesh=al.mesh.Delaunay(),\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])\n", "\n", "fit = al.FitImaging(dataset=dataset, tracer=tracer)\n", "\n", "fit_plotter = aplt.FitImagingPlotter(fit=fit)\n", "fit_plotter.subplot_fit()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Voronoi__\n", "\n", "The pixelization mesh which tests have revealed performs best is the `VoronoiNN` object, which uses a Voronoi\n", "mesh with a technique called natural neighbour interpolation (full details are provided in the **HowToLens**\n", "tutorials).\n", "\n", "I recommend users use this pixelization, however it requires a c library to be installed, thus it is\n", "not the default pixelization used in this tutorial.\n", "\n", "If you want to use this pixelization, checkout the installation instructions here:\n", "\n", "https://github.com/Jammy2211/PyAutoArray/tree/main/autoarray/util/nn\n", "\n", "The code below is commented out because it will not run on your computer, unless you install the c library." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixelization = al.Pixelization(\n", " image_mesh=al.image_mesh.Overlay(shape=(40, 40)),\n", " mesh=al.mesh.Delaunay(),\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])\n", "\n", "fit = al.FitImaging(dataset=dataset, tracer=tracer)\n", "\n", "# Commented out because error will be raised if natural neighbor interpolation is not installed.\n", "# Uncomment if you have installed it!\n", "\n", "# fit_plotter = aplt.FitImagingPlotter(fit=fit)\n", "# fit_plotter.figures_2d_of_planes(plane_index=1, plane_image=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Image Meshes__\n", "\n", "The `Overlay` image mesh above placed a uniform grid of Cartesian coordinates over the image-plane, producing a \n", "source-plane pixelization which adapted to the magnification pattern of the mass-model. The source reconstruction is\n", "good, but could be improved by placing more source-pixels in the source's brightest regions (which typically do not\n", "correspond to the highest magnification regions).\n", "\n", "**PyAutoLens** has alternative image-mesh objects that adapt the source-plane pixelization to the morphology of the \n", "reconstructed unlensed source galaxy. \n", "\n", "This produces the desirable result that high resolution is dedicated to a source's bright, irregular and clumpy \n", "features (e.g. star forming clumps) and fewer pixels to the outskirts which have no measureable signal. \n", "\n", "This type of image-mesh is the recommended approach to source analysis in **PyAutoLens**. It requires the advanced \n", "adaptive pixelization features, which are recommend for experienced PyAutoLens users and \n", "described at `autolens_workspace/*/imaging/advanced/chaining/pix_adapt`.\n", "\n", "__Wrap Up__\n", "\n", "This script has given a brief overview of pixelizations.\n", "\n", "A full descriptions of this feature, including an example of how to use pixelizations in lens modeling, \n", "is given in the `pixelization` example:\n", "\n", "https://github.com/Jammy2211/autolens_workspace/blob/release/notebooks/imaging/modeling/features/pixelization.ipynb\n", "\n", "In chapter 4 of the **HowToLens** lectures we fully cover all aspects of using pixelizations, including:\n", "\n", " - How the source reconstruction determines the flux-values of the source it reconstructs.\n", " - The Bayesian framework employed to choose the appropriate level of regularization and avoid overfitting noise.\n", " - Unphysical lens model solutions that often arise when using pixelizations\n", " - Advanced pixelizations that adapt their properties (e.g. the source pixel locations) to the source galaxy being \n", " reconstructed." ] }, { "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 }