{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Tutorial 3: Lens and Source\n", "===========================\n", "\n", "In this tutorial, we demonstrate search chaining using three searches to fit strong lens `Imaging` which includes the\n", "lens galaxy's light.\n", "\n", "The crucial point to note is that for many lenses the lens galaxy's light can be fitted and subtracted reasonably \n", "well before we attempt to fit the source galaxy. This makes sense, as fitting the lens's light (which is an elliptical\n", "blob of light in the centre of the imaging) looks nothing like the source's light (which is a ring of light)! Formally,\n", "we would say that these two model components (the lens's light and source's light) are not covariate.\n", "\n", "So, as a newly trained lens modeler, what does the lack of covariance between these parameters make you think?\n", "Hopefully, you're thinking, why should I bother fitting the lens and source galaxy simultaneously? Surely we can\n", "find the right regions of non-linear parameter space by fitting each separately first? This is what we're going to do\n", "in this tutorial, using a pipeline composed of a modest 3 searches:\n", "\n", " 1) Fit the lens galaxy's light, ignoring the source.\n", " 2) Fit the source-galaxy's light (and therefore lens galaxy's mass), ignoring the len`s light.\n", " 3) Fit both simultaneously, using these results to initialize our starting location in parameter space.\n", "\n", "Of course, given that we do not care for the errors in searches 1 and 2, we will set up our non-linear search to\n", "perform sampling as fast as possible!\n", "\n", "__Dated Tutorial__\n", "\n", "This example tutorial was written ~4 years ago, when **PyAutoLens** was in its infancy and had a number of limitations:\n", "\n", " - The non-linear search used MultiNest or dynesty, which were less reliable (e.g. more likely to infer a local maxima\n", " for complex lens models) and less efficient than Nautilus.\n", "\n", " - Linear light profiles and techniques like a Multi-Gaussian Expansion were not available.\n", "\n", "With all the new features added to **PyAutoLens** since, we no longer recommend that one breaks down the fitting of\n", "the lens and source galaxy's light into separate searches, as perform in this search chaining example. Instead, we\n", "would recommend you fit the lens and source simultaneously, using linear light profiles to make the model simpler\n", "or a Multi-Gaussian Expansion.\n", "\n", "However, the example is still useful for demonstrating the core concepts of search chaining, which is still vital\n", "for fitting complex lens model. Therefore, we recommend you still read through this tutorial and try to get a good\n", "understanding of how search chaining works, but bear in mind that the example is a little dated and we now recommend\n", "you fit the lens and source simultaneously!" ] }, { "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 autofit as af\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 strong lensing data, where:\n", "\n", " - The lens galaxy's light is an `Sersic`.\n", " - The lens galaxy's total mass distribution is an `Isothermal` and `ExternalShear`.\n", " - The source galaxy's light is an `Exponential`.\n", " \n", "This image was fitted throughout chapter 2." ] }, { "cell_type": "code", "metadata": {}, "source": [ "dataset_name = \"lens_sersic\"\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": [ "__Paths__\n", "\n", "All three searches will use the same `path_prefix`, so we write it here to avoid repetition." ] }, { "cell_type": "code", "metadata": {}, "source": [ "path_prefix = path.join(\"howtolens\", \"chapter_3\", \"tutorial_3_lens_and_source\")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Masking (Search 1)__\n", "\n", "We need to choose our mask for the analysis. Given we are only fitting the lens light we have two options: \n", "\n", " - A circular mask that does not remove the source's light from the fit, assuming the lens light model will still be \n", " sufficiently accurate to reveal the source in the second search.\n", " - An 'anti-annular' mask that removes the source's ring of light.\n", "\n", "In this example, we will use the anti-annular mask to demonstrate that we can change the mask used by each search in a \n", "chain of non-linear searches." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mask = al.Mask2D.circular_anti_annular(\n", " shape_native=dataset.shape_native,\n", " pixel_scales=dataset.pixel_scales,\n", " inner_radius=0.8,\n", " outer_radius=2.2,\n", " outer_radius_2=3.0,\n", ")\n", "\n", "dataset = dataset.apply_mask(mask=mask)\n", "\n", "dataset_plotter = aplt.ImagingPlotter(\n", " dataset=dataset, visuals_2d=aplt.Visuals2D(mask=mask)\n", ")\n", "dataset_plotter.subplot_dataset()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Model + Search + Analysis + Model-Fit (Search 1)__\n", "\n", "Search 1 fits a lens model where:\n", "\n", " - The lens galaxy's light is a parametric linear `Sersic` bulge [6 parameters].\n", " \n", " - The lens galaxy's mass and source galaxy are omitted.\n", "\n", "The number of free parameters and therefore the dimensionality of non-linear parameter space is N=6.\n", "\n", "__Notes__\n", "\n", "We use linear light profiles througout this script, given that the model is quite complex and this helps\n", "simplify it." ] }, { "cell_type": "code", "metadata": {}, "source": [ "model_1 = af.Collection(\n", " galaxies=af.Collection(lens=af.Model(al.Galaxy, redshift=0.5, bulge=al.lp.Sersic)),\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `info` attribute shows the model in a readable format." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(model_1.info)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Search + Analysis + Model-Fit (Search 1)__" ] }, { "cell_type": "code", "metadata": {}, "source": [ "analysis_1 = al.AnalysisImaging(dataset=dataset)\n", "\n", "search_1 = af.Nautilus(\n", " path_prefix=path_prefix,\n", " name=\"search[1]_light[bulge]\",\n", " unique_tag=dataset_name,\n", " n_live=75,\n", " number_of_cores=4,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Run Time__\n", "\n", "It is good practise to always check the `log_likelihood_function` run time before starting the non-linear search. \n", "It will be similar to the value we saw in the previous chapter." ] }, { "cell_type": "code", "metadata": {}, "source": [ "run_time_dict, info_dict = analysis_1.profile_log_likelihood_function(\n", " instance=model_1.random_instance()\n", ")\n", "\n", "print(f\"Log Likelihood Evaluation Time (second) = {run_time_dict['fit_time']}\")\n", "print(\n", " \"Estimated Run Time Upper Limit (seconds) = \",\n", " (run_time_dict[\"fit_time\"] * model_1.total_free_parameters * 10000)\n", " / search_1.number_of_cores,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the search." ] }, { "cell_type": "code", "metadata": {}, "source": [ "\n", "result_1 = search_1.fit(model=model_1, analysis=analysis_1)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Result (Search 1)__\n", "\n", "The results which are used for prior passing are summarized in the `info` attribute." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(result_1.info)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Masking (Search 2)__\n", "\n", "Search 2 we are only fitting the source's light, thus we can apply an annular mask that removes regions of the\n", "image that contained only the lens's 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.6,\n", " outer_radius=2.4,\n", ")\n", "\n", "dataset = dataset.apply_mask(mask=mask)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Model + Search + Analysis + Model-Fit (Search 2)__\n", "\n", "Search 2 fits a lens model where:\n", "\n", " - The lens galaxy's light is a linear `Sersic` bulge [Parameters fixed to results of search 1].\n", " \n", " - The lens galaxy's total mass distribution is an `Isothermal` with `ExternalShear` [7 parameters].\n", " \n", " - The source galaxy's light is a parametric linear `Sersic` [6 parameters].\n", "\n", "The number of free parameters and therefore the dimensionality of non-linear parameter space is N=13.\n", "\n", "Search 2, we fit the source-`galaxy's light and fix the lens light model to the model inferred in search 1, \n", "ensuring the image we has the foreground lens subtracted. We do this below by passing the lens light as an `instance` \n", "object.\n", "\n", "By passing an `instance`, we are telling **PyAutoLens** that we want it to pass the maximum log likelihood result of \n", "that search and use those parameters as fixed values in the model. The model parameters passed as an `instance` are not \n", "free parameters fitted for by the non-linear search, thus this reduces the dimensionality of the non-linear search \n", "making model-fitting faster and more reliable. \n", " \n", "Thus, search 2 includes the lens light model from search 1, but it is completely fixed during the model-fit!\n", "\n", "We also use the centre of the `bulge` to initialize the priors on the lens's `mass`." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mass = af.Model(al.mp.Isothermal)\n", "mass.centre_0 = result_1.model.galaxies.lens.bulge.centre_0\n", "mass.centre_1 = result_1.model.galaxies.lens.bulge.centre_1\n", "\n", "model_2 = af.Collection(\n", " galaxies=af.Collection(\n", " lens=af.Model(\n", " al.Galaxy,\n", " redshift=0.5,\n", " bulge=result_1.instance.galaxies.lens.bulge,\n", " mass=mass,\n", " shear=al.mp.ExternalShear,\n", " ),\n", " source=af.Model(al.Galaxy, redshift=1.0, bulge=al.lp.Sersic),\n", " ),\n", ")\n", "\n", "analysis_2 = al.AnalysisImaging(dataset=dataset)\n", "\n", "search_2 = af.Nautilus(\n", " path_prefix=path_prefix,\n", " name=\"search[2]_mass[sie]_source[bulge]\",\n", " unique_tag=dataset_name,\n", " n_live=100,\n", " number_of_cores=4,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Run Time__\n", "\n", "The run-time of the fit should be noticeably faster than the previous search, but because the smaller mask means the\n", "likelihood function is evaluated faster and because prior passing ensures the search samples parameter space faster." ] }, { "cell_type": "code", "metadata": {}, "source": [ "run_time_dict, info_dict = analysis_2.profile_log_likelihood_function(\n", " instance=model_2.random_instance()\n", ")\n", "\n", "print(f\"Log Likelihood Evaluation Time (second) = {run_time_dict['fit_time']}\")\n", "print(\n", " \"Estimated Run Time Upper Limit (seconds) = \",\n", " (run_time_dict[\"fit_time\"] * model_2.total_free_parameters * 10000)\n", " / search_2.number_of_cores,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the search." ] }, { "cell_type": "code", "metadata": {}, "source": [ "result_2 = search_2.fit(model=model_2, analysis=analysis_2)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Result (Search 2)__\n", "\n", "The results which are used for prior passing are summarized in the `info` attribute." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(result_2.info)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Masking (Search 3)__\n", "\n", "Search 3 we fit the lens and source, therefore we will use a large circular mask." ] }, { "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", "dataset = dataset.apply_mask(mask=mask)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Model + Search + Analysis + Model-Fit (Search 3)__\n", "\n", "Search 3 fits a lens model where:\n", "\n", " - The lens galaxy's light is a linear `Sersic` bulge [6 Parameters: priors initialized from search 1].\n", " \n", " - The lens galaxy's total mass distribution is an `Isothermal` with `ExternalShear` [7 parameters: priors\n", " initialized from search 2].\n", " \n", " - The source galaxy's light is a parametric linear `Sersic` [6 parameters: priors initialized from search 2].\n", "\n", "The number of free parameters and therefore the dimensionality of non-linear parameter space is N=23.\n", "\n", "There isn't a huge amount to say about this search, we have initialized the priors on all of our models parameters\n", "and the only thing that is left to do is fit for all model components simultaneously, with slower Nautilus settings\n", "that will give us more accurate parameter values and errors." ] }, { "cell_type": "code", "metadata": {}, "source": [ "model_3 = af.Collection(\n", " galaxies=af.Collection(\n", " lens=af.Model(\n", " al.Galaxy,\n", " redshift=0.5,\n", " bulge=result_1.model.galaxies.lens.bulge,\n", " mass=result_2.model.galaxies.lens.mass,\n", " ),\n", " source=af.Model(\n", " al.Galaxy, redshift=1.0, bulge=result_2.model.galaxies.source.bulge\n", " ),\n", " ),\n", ")\n", "\n", "analysis_3 = al.AnalysisImaging(dataset=dataset)\n", "\n", "search_3 = af.Nautilus(\n", " path_prefix=path_prefix,\n", " name=\"search[3]_light[bulge]_mass[sie]_source[bulge]\",\n", " unique_tag=dataset_name,\n", " n_live=150,\n", " number_of_cores=4,\n", ")\n", "\n", "result_3 = search_3.fit(model=model_3, analysis=analysis_3)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Result (Search 3)__\n", "\n", "The final results are summarized in the `info` attribute." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(result_3.info)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Wrap Up__\n", "\n", "And there we have it, a sequence of searches that breaks modeling the lens and source galaxy into 3 simple searches. \n", "This approach is much faster than fitting the lens and source simultaneously from the beginning. Instead of asking you \n", "questions at the end of this chapter`s tutorials, I'm going to ask questions which I then answer. This will hopefully \n", "get you thinking about how to approach pipeline writing.\n", "\n", " 1) Can this pipeline really be generalized to any lens? Surely the radii of the masks depends on the lens and source \n", " galaxies?\n", "\n", "Whilst this is true, we chose mask radii above that are `excessive` and masks out a lot more of the image than just \n", "the source (which, in terms of run-time, is desirable). Thus, provided you know the Einstein radius distribution of \n", "your lens sample, you can choose mask radii that will masks out every source in your sample adequately (and even if \n", "some of the source is still there, who cares? The fit to the lens galaxy will be okay).\n", "\n", "However, the template pipelines provided on the `autolens_workspace` simply use circular masks for every search and do\n", "not attempt to use different masks for the lens light fit and source fit. This is to keep things simple (at the expense\n", "of slower run times). It is up to you if you want to adapt these scripts to try and use more specific masking strategies.\n", "\n", "__Dated Tutorial__\n", "\n", "In fact, we now strongly recommend that you do not change masks between each search when using search chaining. \n", "This is because it is very fiddly, and can waste a lot of your time refining masks to ensure they are suitable for\n", "each lens. We recommend you always just use a large circular mask which is big enough to include the entire lens and \n", "source of all lenses in your sample. This will save you a lot of time and means lens modeling can be automated much\n", "easier.\n", "\n", "Building on the discussion above, a known limitation of using a pipeline which fits the lens light first, then the\n", "source, is that it will do a poor job deblending the lens and source light if the Einstein radius is low. This often\n", "leads the mass model to infer incorrect solutions which fit residuals from the lens light subtraction.\n", "\n", "This is why, given all the improvements to autolens, we now recommend that you do not use this pipeline and instead\n", "always begin by fitting the lens and source simultaneously. This can use linear light profiles of a Multi-Gaussian\n", "Expansion. " ] }, { "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 }