{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Tutorial 11: Adaptive Regularization\n", "====================================\n", "\n", "In tutorial 7, we discussed why the `Constant` regularization scheme was sub-optimal. Different regions of the source\n", "demand different levels of regularization, motivating a regularization scheme which adapts to the reconstructed\n", "source's surface brightness.\n", "\n", "This raises the same question as before, how do we adapt our regularization scheme to the source before we've\n", "reconstructed it? Just like in the last tutorial, we'll use a model image of a strongly lensed source from a previous\n", "model that we've begun calling the `adapt-image`." ] }, { "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", ")\n", "\n", "mask = al.Mask2D.circular(\n", " shape_native=dataset.shape_native,\n", " pixel_scales=dataset.pixel_scales,\n", " sub_size=2,\n", " radius=3.0,\n", ")\n", "\n", "dataset = dataset.apply_mask(mask=mask)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Convenience Function__\n", "\n", "We are going to fit the image using a magnification based grid. \n", "\n", "To perform the fits, we'll use a convenience function to fit the lens data we simulated above." ] }, { "cell_type": "code", "metadata": {}, "source": [ "\n", "\n", "def fit_via_source_galaxy_from(dataset, source_galaxy, adapt_images=None):\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.from_galaxies(galaxies=[lens_galaxy, source_galaxy])\n", "\n", " return al.FitImaging(dataset=dataset, tracer=tracer, adapt_images=adapt_images)\n" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the magnification based source to fit this data." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixelization = al.Pixelization(\n", " image_mesh=al.image_mesh.Overlay(shape=(30, 30)),\n", " mesh=al.mesh.Delaunay(),\n", " regularization=al.reg.Constant(coefficient=3.3),\n", ")\n", "\n", "source_magnification = al.Galaxy(redshift=1.0, pixelization=pixelization)\n", "\n", "fit = fit_via_source_galaxy_from(dataset=dataset, source_galaxy=source_magnification)\n", "\n", "include = aplt.Include2D(\n", " mask=True, mapper_image_plane_mesh_grid=True, mapper_source_plane_mesh_grid=True\n", ")\n", "\n", "fit_plotter = aplt.FitImagingPlotter(fit=fit, include_2d=include)\n", "fit_plotter.subplot_fit()\n", "fit_plotter.figures_2d_of_planes(plane_index=1, plane_image=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fit looks just like it did in the previous tutorials (residuals in the centre due to a lack of source pixels). \n", "\n", "Lets quickly remind ourselves that the effective regularization weight of each source pixel is our input coefficient \n", "value of 3.3." ] }, { "cell_type": "code", "metadata": {}, "source": [ "inversion_plotter = fit_plotter.inversion_plotter_of_plane(plane_index=1)\n", "inversion_plotter.figures_2d_of_pixelization(\n", " pixelization_index=0, regularization_weights=True\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Adaptive Regularization__\n", "\n", "Lets now look at adaptive regularization in action, by setting up a adapt-image and using the `AdaptiveBrightness` \n", "regularization scheme. \n", "\n", "This introduces additional parameters, that are explained below." ] }, { "cell_type": "code", "metadata": {}, "source": [ "adapt_image = fit.model_image.binned.slim\n", "\n", "pixelization = al.Pixelization(\n", " image_mesh=al.image_mesh.Overlay(shape=(30, 30)),\n", " mesh=al.mesh.Delaunay(),\n", " regularization=al.reg.AdaptiveBrightness(\n", " inner_coefficient=0.005, outer_coefficient=1.9, signal_scale=3.0\n", " ),\n", ")\n", "\n", "source_adaptive_regularization = al.Galaxy(\n", " redshift=1.0,\n", " pixelization=pixelization,\n", " adapt_galaxy_image=adapt_image,\n", ")\n", "\n", "adapt_images = al.AdaptImages(\n", " galaxy_image_dict={source_adaptive_regularization: adapt_image}\n", ")\n", "\n", "fit = fit_via_source_galaxy_from(\n", " dataset=dataset,\n", " source_galaxy=source_adaptive_regularization,\n", " adapt_images=adapt_images,\n", ")\n", "\n", "inversion_plotter = fit_plotter.inversion_plotter_of_plane(plane_index=1)\n", "inversion_plotter.figures_2d_of_pixelization(\n", " pixelization_index=0, reconstruction=True, regularization_weights=True\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, we now have a variable regularization scheme. \n", "\n", "The regularization of the source's brightest regions is much lower than that of its outer regions. As discussed \n", "before, this is what we want. Lets quickly check that this does, indeed, increase the Bayesian log evidence:" ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(\"Evidence using constant regularization. \", 4216)\n", "print(\"Evidence using adaptive regularization. \", fit.log_evidence)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yes, it does! \n", "\n", "Combining the adaptive mesh and regularization will only further benefit lens modeling!\n", "\n", "However, as shown below, we don't fit the source as well as the morphology based mesh did in the last chapter. \n", "This is because although the adaptive regularization scheme improves the fit, the magnification based \n", "mesh simply does not have sufficient resolution to resolve the source's cuspy central light." ] }, { "cell_type": "code", "metadata": {}, "source": [ "fit_plotter = aplt.FitImagingPlotter(fit=fit, include_2d=include)\n", "fit_plotter.subplot_fit()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__How does adaptive regularization work?__\n", "\n", "For every source-pixel, we have a mapping between that pixel and a set of pixels in the adapt-image. Therefore, for \n", "every source-pixel, if we sum the values of all adapt-image pixels that map to it we get an estimate of how much of \n", "the lensed source's signal we expect will be reconstructed. We call this each pixel's `pixel signal`.\n", "\n", "If a source-pixel has a higher pixel-signal, we anticipate that it`ll reconstruct more flux and we use this information \n", "to regularize it less. Conversely, if the pixel-signal is close to zero, the source pixel will reconstruct near-zero \n", "flux and regularization will smooth over these pixels by using a high regularization coefficient.\n", "\n", "This works as follows:\n", "\n", " 1) For every source pixel, compute its pixel-signal, the summed flux of all corresponding image-pixels in the \n", " adapt-galaxy-image.\n", " \n", " 2) Divide every pixel-signal by the number of image-pixels that map directly to that source-pixel. In doing so, all \n", " pixel-signals are 'relative'. This means that source-pixels which by chance map to more image-pixels than their \n", " neighbors will not have a higher pixel-signal, and visa versa. This ensures the specific pixelization\n", " does impact the adaptive regularization pattern.\n", " \n", " 3) Divide the pixel-signals by the maximum pixel signal so that they range between 0.0 and 1.0.\n", " \n", " 4) Raise these values to the power of the adapt-parameter `signal_scale`. For a `signal_scale` of 0.0, all \n", " pixels will therefore have the same final pixel-scale. As the `signal_scale` increases, a sharper transition of \n", " pixel-signal values arises between regions with high and low pixel-signals.\n", " \n", " 5) Compute every source pixel's effective regularization coefficient as:\n", " \n", " (inner_coefficient * pixel_signals + outer_coefficient * (1.0 - pixel_signals)) ** 2.0\n", " \n", " This uses two regularization coefficient parameters, one which is applied to pixels with high pixel-signals and one to \n", " pixels with low pixel-signals. Thus, pixels in the inner regions of the source may be given a lower level of \n", " regularization than pixels further away, as desired.\n", "\n", "Thus, we now adapt our regularization scheme to the source's surface brightness. Where its brighter (and therefore \n", "has a steeper flux gradient) we apply a lower level of regularization than further out. Furthermore, in the edges of \n", "the source-plane where no source-flux is present we will assume a high regularization coefficient that smooths over \n", "all of the source-pixels.\n", "\n", "Try looking at a couple of extra solutions which use with different inner and outer regularization coefficients or \n", "signal scales. I doubt you'll notice a lot change visually, but the log evidence certainly has a lot of room for \n", "manuveur with different values.\n", "\n", "You may find solutions that raise an `InversionException`. These solutions mean that the matrix used during the \n", "linear algebra calculation was ill-defined, and could not be inverted. These solutions are removed \n", "during lens modeling." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixelization = al.Pixelization(\n", " image_mesh=al.image_mesh.Overlay(shape=(30, 30)),\n", " mesh=al.mesh.Delaunay(),\n", " regularization=al.reg.AdaptiveBrightness(\n", " inner_coefficient=0.001, outer_coefficient=0.2, signal_scale=2.0\n", " ),\n", ")\n", "\n", "source_adaptive_regularization = al.Galaxy(\n", " redshift=1.0,\n", " pixelization=pixelization,\n", " adapt_galaxy_image=adapt_image,\n", ")\n", "\n", "adapt_images = al.AdaptImages(\n", " galaxy_image_dict={source_adaptive_regularization: adapt_image}\n", ")\n", "\n", "fit = fit_via_source_galaxy_from(\n", " dataset=dataset,\n", " source_galaxy=source_adaptive_regularization,\n", " adapt_images=adapt_images,\n", ")\n", "\n", "fit_plotter = aplt.FitImagingPlotter(fit=fit, include_2d=include)\n", "fit_plotter.subplot_fit()\n", "\n", "inversion_plotter = fit_plotter.inversion_plotter_of_plane(plane_index=1)\n", "inversion_plotter.figures_2d_of_pixelization(\n", " pixelization_index=0, reconstruction=True, regularization_weights=True\n", ")\n", "\n", "print(\"Evidence using adaptive regularization. \", fit.log_evidence)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Wrap Up__\n", "\n", "To end, lets consider what this adaptive regularization scheme means in the context of maximizing the Bayesian\n", "evidence. In the previous tutorial, we noted that by using a brightness-based adaptive pixelization we increased \n", "the Bayesian evidence by allowing for new solutions which fit the data user fewer source pixels; the key criteria \n", "in making a source reconstruction 'more simple' and 'less complex'.\n", "\n", "As you might of guessed, adaptive regularization increases the Bayesian log evidence by making the source \n", "reconstruction simpler:\n", "\n", " 1) Reducing regularization in the source's brightest regions produces a `simpler` solution in that we are not \n", " over-smoothing our reconstruction of its brightest regions.\n", " \n", " 2) Increasing regularization in the outskirts produces a simpler solution by correlating more source-pixels, \n", " effectively reducing the number of pixels used by the reconstruction.\n", "\n", "Together, brightness based pixelization's and regularization allow us to find the objectively `simplest` source \n", "solution possible and therefore ensure that our Bayesian evidence has a well defined maximum value. This was not the \n", "case for magnification based pixelization's and constant regularization schemes." ] }, { "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 }