{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Overview: Multi-Wavelength\n", "--------------------------\n", "\n", "**PyAutoLens** supports the analysis of multiple datasets simultaneously, including many CCD imaging datasets\n", "observed at different wavebands (e.g. red, blue, green) and combining imaging and interferometer datasets.\n", "\n", "This enables multi-wavelength lens modeling, where the color of the lens and source galaxies vary across the datasets.\n", "\n", "Multi-wavelength lens modeling offers a number of advantages:\n", "\n", "- It provides a wealth of additional information to fit the lens model, especially if the source changes its\n", "appearance across wavelength.\n", "\n", "- It overcomes challenges associated with the lens and source galaxy emission blending with one another, as their\n", " brightness depends differently on wavelength.\n", "\n", "- Instrument systematic effects, for example an uncertain PSF, will impact the model less because they vary across\n", " each dataset." ] }, { "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 autofit as af\n", "import autolens as al\n", "import autolens.plot as aplt\n", "from os import path\n", "import numpy as np" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Colors__\n", "\n", "For multi-wavelength imaging datasets, we begin by defining the colors of the multi-wavelength images. \n", "\n", "For this overview we use only two colors, green (g-band) and red (r-band), but extending this to more datasets\n", "is straight forward." ] }, { "cell_type": "code", "metadata": {}, "source": [ "color_list = [\"g\", \"r\"]" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Pixel Scales__\n", "\n", "Every dataset in our multi-wavelength observations can have its own unique pixel-scale." ] }, { "cell_type": "code", "metadata": {}, "source": [ "pixel_scales_list = [0.08, 0.12]" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Dataset__\n", "\n", "Multi-wavelength imaging datasets do not use any new objects or class in **PyAutoLens**.\n", "\n", "We simply use lists of the classes we are now familiar with, for example the `Imaging` class." ] }, { "cell_type": "code", "metadata": {}, "source": [ "dataset_type = \"multi\"\n", "dataset_label = \"imaging\"\n", "dataset_name = \"lens_sersic\"\n", "\n", "dataset_path = path.join(\"dataset\", dataset_type, dataset_label, dataset_name)\n", "\n", "dataset_list = [\n", " al.Imaging.from_fits(\n", " data_path=path.join(dataset_path, f\"{color}_data.fits\"),\n", " psf_path=path.join(dataset_path, f\"{color}_psf.fits\"),\n", " noise_map_path=path.join(dataset_path, f\"{color}_noise_map.fits\"),\n", " pixel_scales=pixel_scales,\n", " )\n", " for color, pixel_scales in zip(color_list, pixel_scales_list)\n", "]" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is what our r-band and g-band observations of this lens system looks like.\n", "\n", "In the r-band, the lens outshines the source, whereas in the g-band the source galaxy is more visible. \n", "\n", "The different variation of the colors of the lens and source across wavelength is a powerful tool for lens modeling,\n", "as it helps deblend the two objects." ] }, { "cell_type": "code", "metadata": {}, "source": [ "for dataset in dataset_list:\n", " dataset_plotter = aplt.ImagingPlotter(dataset=dataset)\n", " dataset_plotter.subplot_dataset()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Mask__\n", "\n", "Define a 3.0\" circular mask, which includes the emission of the lens and source galaxies.\n", "\n", "For multi-wavelength lens modeling, we use the same mask for every dataset whenever possible. This is not absolutely \n", "necessary, but provides a more reliable analysis." ] }, { "cell_type": "code", "metadata": {}, "source": [ "mask_list = [\n", " al.Mask2D.circular(\n", " shape_native=dataset.shape_native, pixel_scales=dataset.pixel_scales, radius=3.0\n", " )\n", " for dataset in dataset_list\n", "]\n", "\n", "\n", "dataset_list = [\n", " dataset.apply_mask(mask=mask) for dataset, mask in zip(dataset_list, mask_list)\n", "]\n", "\n", "for dataset in dataset_list:\n", " dataset_plotter = aplt.ImagingPlotter(dataset=dataset)\n", " dataset_plotter.subplot_dataset()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Analysis__\n", "\n", "We create a list of `AnalysisImaging` objects for every dataset." ] }, { "cell_type": "code", "metadata": {}, "source": [ "analysis_list = [al.AnalysisImaging(dataset=dataset) for dataset in dataset_list]" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now introduce the key new aspect to the **PyAutoLens** multi-dataset API, which is critical to fitting multiple \n", "datasets simultaneously.\n", "\n", "We sum the list of analysis objects to create an overall `CombinedAnalysis` object, which we can use to fit the \n", "multi-wavelength imaging dataset, where:\n", "\n", " - The log likelihood function of this summed analysis class is the sum of the log likelihood functions of each \n", " individual analysis objects (e.g. the fit to each separate waveband).\n", "\n", " - The summing process ensures that tasks such as outputting results to hard-disk, visualization, etc use a \n", " structure that separates each analysis and therefore each dataset." ] }, { "cell_type": "code", "metadata": {}, "source": [ "analysis = sum(analysis_list)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Model__\n", "\n", "We compose an initial lens model as per usual." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# Lens:\n", "\n", "bulge = af.Model(al.lp.Sersic)\n", "\n", "mass = af.Model(al.mp.Isothermal)\n", "mass.centre = (0.0, 0.0)\n", "\n", "shear = af.Model(al.mp.ExternalShear)\n", "\n", "lens = af.Model(\n", " al.Galaxy,\n", " redshift=0.5,\n", " bulge=bulge,\n", " mass=mass,\n", " shear=shear,\n", ")\n", "\n", "# Source:\n", "\n", "bulge = af.Model(al.lp.Sersic)\n", "\n", "source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge)\n", "\n", "# Overall Lens Model:\n", "\n", "model = af.Collection(galaxies=af.Collection(lens=lens, source=source))" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, there is a problem for multi-wavelength datasets. Should the light profiles of the lens's bulge and\n", "source's bulge have the same parameters for each wavelength image?\n", "\n", "The answer is no. At different wavelengths, different stars appear brighter or fainter, meaning that the overall\n", "appearance of the lens and source galaxies will change. \n", "\n", "We therefore allow specific light profile parameters to vary across wavelength and act as additional free\n", "parameters in the fit to each image. \n", "\n", "We do this using the combined analysis object as follows:" ] }, { "cell_type": "code", "metadata": {}, "source": [ "analysis = analysis.with_free_parameters(\n", " model.galaxies.lens.bulge.intensity, model.galaxies.source.bulge.intensity\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this simple overview, this has added two additional free parameters to the model whereby:\n", "\n", " - The lens bulge's intensity is different in both multi-wavelength images.\n", " - The source bulge's intensity is different in both multi-wavelength images.\n", " \n", "It is entirely plausible that more parameters should be free to vary across wavelength (e.g. the lens and source\n", "galaxies `effective_radius` or `sersic_index` parameters). \n", "\n", "This choice ultimately depends on the quality of data being fitted and intended science goal. Regardless, it is clear\n", "how the above API can be extended to add any number of additional free parameters.\n", "\n", "__Search + Model Fit__\n", "\n", "Fitting the model uses the same API we introduced in previous overviews." ] }, { "cell_type": "code", "metadata": {}, "source": [ "search = af.Nautilus(\n", " path_prefix=\"overview\",\n", " name=\"multi_wavelength\",\n", " n_live=200,\n", ")" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result object returned by this model-fit is a list of `Result` objects, because we used a combined analysis.\n", "Each result corresponds to each analysis created above and is there the fit to each dataset at each wavelength." ] }, { "cell_type": "code", "metadata": {}, "source": [ "result_list = search.fit(model=model, analysis=analysis)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting each result's tracer shows that the lens and source galaxies appear different in each result, owning to their \n", "different intensities." ] }, { "cell_type": "code", "metadata": {}, "source": [ "for result in result_list:\n", " tracer_plotter = aplt.TracerPlotter(\n", " tracer=result.max_log_likelihood_tracer, grid=result.grid\n", " )\n", " tracer_plotter.subplot_tracer()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Subplots of each fit show that a good overall fit is achieved to each dataset." ] }, { "cell_type": "code", "metadata": {}, "source": [ "for result in result_list:\n", " fit_plotter = aplt.FitImagingPlotter(\n", " fit=result.max_log_likelihood_fit,\n", " )\n", " fit_plotter.subplot_fit()" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Wavelength Dependence__\n", "\n", "In the example above, a free `intensity` parameter is created for every multi-wavelength dataset. This would add 5+ \n", "free parameters to the model if we had 5+ datasets, quickly making a complex model parameterization.\n", "\n", "We can instead parameterize the intensity of the lens and source galaxies as a user defined function of \n", "wavelength, for example following a relation `y = (m * x) + c` -> `intensity = (m * wavelength) + c`.\n", "\n", "By using a linear relation `y = mx + c` the free parameters are `m` and `c`, which does not scale with the number\n", "of datasets. For datasets with multi-wavelength images (e.g. 5 or more) this allows us to parameterize the variation \n", "of parameters across the datasets in a way that does not lead to a very complex parameter space.\n", "\n", "Below, we show how one would do this for the `intensity` of a lens galaxy's bulge, give three wavelengths corresponding\n", "to a dataset observed in the g, r and I bands." ] }, { "cell_type": "code", "metadata": {}, "source": [ "wavelength_list = [464, 658, 806]\n", "\n", "lens_m = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1)\n", "lens_c = af.UniformPrior(lower_limit=-10.0, upper_limit=10.0)\n", "\n", "source_m = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1)\n", "source_c = af.UniformPrior(lower_limit=-10.0, upper_limit=10.0)\n", "\n", "analysis_list = []\n", "\n", "for wavelength, dataset in zip(wavelength_list, dataset_list):\n", " lens_intensity = (wavelength * lens_m) + lens_c\n", " source_intensity = (wavelength * source_m) + source_c\n", "\n", " analysis_list.append(\n", " al.AnalysisImaging(dataset=dataset).with_model(\n", " model.replacing(\n", " {\n", " model.galaxies.lens.bulge.intensity: lens_intensity,\n", " model.galaxies.source.bulge.intensity: source_intensity,\n", " }\n", " )\n", " )\n", " )" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Same Wavelength Datasets__\n", "\n", "The above API can fit multiple datasets which are observed at the same wavelength.\n", "\n", "For example, this allows the analysis of images of a galaxy before they are combined to a single frame via the \n", "multidrizzling data reduction process to remove correlated noise in the data.\n", "\n", "The pointing of each observation, and therefore centering of each dataset, may vary in an unknown way. This\n", "can be folded into the model and fitted for as follows:\n", "\n", "TODO : ADD CODE EXAMPLE.\n", "\n", "__Interferometry and Imaging__\n", "\n", "The above API can combine modeling of imaging and interferometer datasets \n", "(see `autolens_workspace/*/multi/modeling/imaging_and_interferometer.ipynb` for an example script showing \n", "this in full).\n", "\n", "Below are mock strong lens images of a system observed at a green wavelength (g-band) and with an interferometer at\n", "sub millimeter wavelengths. \n", "\n", "A number of benefits are apparent if we combine the analysis of both datasets at both wavelengths:\n", "\n", " - The lens galaxy is invisible at sub-mm wavelengths, making it straight-forward to infer a lens mass model by\n", " fitting the source at submm wavelengths.\n", " \n", " - The source galaxy appears completely different in the g-band and at sub-millimeter wavelengths, providing a lot\n", " more information with which to constrain the lens galaxy mass model.\n", " \n", "__Linear Light Profiles__\n", "\n", "The modeling overview example described linear light profiles, where the `intensity` parameters of all parametric \n", "components are solved via linear algebra every time the model is fitted using a process called an inversion. \n", "\n", "These profiles are particular powerful when combined with multi-wavelength datasets, because the linear algebra\n", "will solve for the `intensity` value in each individual dataset separately. \n", "\n", "This means that the `intensity` value of all of the galaxy light profiles in the model vary across the multi-wavelength\n", "datasets, but the dimensionality of the model does not increase as more datasets are fitted.\n", "\n", "A full example is given in the `linear_light_profiles` example:\n", "\n", "https://github.com/Jammy2211/autolens_workspace/blob/release/notebooks/multi/modeling/features/linear_light_profiles.ipynb\n", "\n", "__Wrap Up__\n", "\n", "The `multi `_ package\n", "of the `autolens_workspace `_ contains numerous example scripts for performing\n", "multi-wavelength modeling and simulating strong lenses with multiple datasets." ] }, { "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 }