PyAutoLens is software for analysing strong gravitational lenses, an astrophysical phenomenon where a galaxy appears multiple times because its light is bent by the gravitational field of an intervening foreground lens galaxy.
Here is a schematic of a strong gravitational lens:
Credit: F. Courbin, S. G. Djorgovski, G. Meylan, et al., Caltech / EPFL / WMKO https://www.astro.caltech.edu/~george/qsolens/
This notebook gives an overview of PyAutoLens's features and API.
Lets first import autolens, its plotting module and the other libraries we'll need.
%matplotlib inline import autolens as al import autolens.plot as aplt import matplotlib.pyplot as plt from os import path
First, we illustrate lensing calculations in PyAutoLens by creating an an image of a strong lens.
To describe the deflection of light, PyAutoLens uses
Grid2D data structures, which are two-dimensional
Cartesian grids of (y,x) coordinates.
Below, we make and plot a uniform Cartesian grid:
grid_2d = al.Grid2D.uniform( shape_native=(150, 150), pixel_scales=0.05, # <- The pixel-scale describes the conversion from pixel units to arc-seconds. ) grid_2d_plotter = aplt.Grid2DPlotter(grid=grid_2d) grid_2d_plotter.figure_2d()
Our aim is to create an image of the source galaxy after its light has been deflected by the mass of the foreground
lens galaxy. We therefore need to ray-trace the
Grid2D's coordinates from the 'image-plane' to the 'source-plane'.
We therefore need analytic functions representing a galaxy's light and mass distributions. For this, PyAutoLens
Profile objects, for example an ellipical sersic
LightProfile object which represents a light distribution:
sersic_light_profile = al.lp.Sersic( centre=(0.0, 0.0), ell_comps=(0.2, 0.1), intensity=0.005, effective_radius=2.0, sersic_index=4.0, )
By passing this profile the
Grid2D, we can evaluate the light emitted at every (y,x) coordinate on the
create an image of the
image_2d = sersic_light_profile.image_2d_from(grid=grid_2d) plt.imshow(image_2d.native) # The use of 'native' is described at the start of the HowToLens tutorials.
<matplotlib.image.AxesImage at 0x7f38fed1f520>
The PyAutoLens plot module provides methods for plotting objects and their properties, like the image of
Note how, unlike the matplotlib method above, this figure is displayed with axis units of arc-seconds, a colorbar, labels, a title, etc. The PyAutoLens plot module takes care of all the heavy lifting that comes with making figures!
light_profile_plotter = aplt.LightProfilePlotter( light_profile=sersic_light_profile, grid=grid_2d ) light_profile_plotter.figures_2d(image=True)
MassProfile objects to represent a galaxy's mass distribution and perform ray-tracing
Below we create an elliptical isothermal
MassProfile and compute its deflection angles on our Cartesian grid, where
the deflection angles describe how the lens galaxy's mass bends the source's light:
isothermal_mass_profile = al.mp.Isothermal( centre=(0.0, 0.0), ell_comps=(0.1, 0.0), einstein_radius=1.6 ) deflections = isothermal_mass_profile.deflections_yx_2d_from(grid=grid_2d)
The deflection angles are easily plotted using the PyAutoLens plot module.
mass_profile_plotter = aplt.MassProfilePlotter( mass_profile=isothermal_mass_profile, grid=grid_2d ) mass_profile_plotter.figures_2d( deflections_y=True, deflections_x=True )
Many other lensing quantities are easily plotted with the
For readers unsure what the plotted quantities below mean, they are described in detail in chapter 1 of the HowToLens Jupyter notebook lectures:
mass_profile_plotter.figures_2d( convergence=True, magnification=True, )
In PyAutoLens, a
Galaxy object is a collection of
MassProfile objects at a given redshift.
The code below creates two galaxies representing the lens and source galaxies shown in the strong lensing diagram above.
lens_galaxy = al.Galaxy( redshift=0.5, light=sersic_light_profile, mass=isothermal_mass_profile ) source_light_profile = al.lp.Exponential( centre=(0.3, 0.2), ell_comps=(0.1, 0.0), intensity=0.1, effective_radius=0.5 ) source_galaxy = al.Galaxy(redshift=1.0, light=source_light_profile)
We can use a
GalaxyPlotter to plot the properties of the lens and source galaxies.
lens_galaxy_plotter = aplt.GalaxyPlotter(galaxy=lens_galaxy, grid=grid_2d) lens_galaxy_plotter.figures_2d(image=True, deflections_y=True, deflections_x=True) source_galaxy_plotter = aplt.GalaxyPlotter(galaxy=source_galaxy, grid=grid_2d) source_galaxy_plotter.figures_2d(image=True)
The geometry of the strong lens system depends on the cosmological distances between the Earth, the lens galaxy and
the source galaxy. It therefore depends on the redshifts of the
By passing these
Galaxy objects to a
Tracer, PyAutoLens uses these galaxy redshifts and a cosmological model
to create the appropriate strong lens system.
tracer = al.Tracer.from_galaxies( galaxies=[lens_galaxy, source_galaxy], cosmology=al.cosmo.Planck15() ) image_2d = tracer.image_2d_from(grid=grid_2d) tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid_2d) tracer_plotter.figures_2d(image=True)
TracerPlotter includes the
MassProfile quantities we plotted previously, which can be plotted as a subplot
that plots all these quantities simultaneously.
All of the objects introduced above are extensible.
Galaxy objects can take many
If the galaxies are at different redshifts a strong lensing system with multiple lens planes will be created, performing complex multi-plane ray-tracing calculations.
To finish, lets create a
Tracer with 3 galaxies at 3 different redshifts, forming a system with two distinct Einstein
rings! The mass distribution of the first galaxy has separate components for its stellar mass and dark matter, where
the stellar components use a
LightAndMassProfile via the
lens_galaxy_0 = al.Galaxy( redshift=0.5, bulge=al.lmp.Sersic( centre=(0.0, 0.0), ell_comps=(0.0, 0.05), intensity=0.5, effective_radius=0.3, sersic_index=3.5, mass_to_light_ratio=0.6, ), disk=al.lmp.Exponential( centre=(0.0, 0.0), ell_comps=(0.0, 0.1), intensity=1.0, effective_radius=2.0, mass_to_light_ratio=0.2, ), dark=al.mp.NFWSph(centre=(0.0, 0.0), kappa_s=0.08, scale_radius=30.0), ) lens_galaxy_1 = al.Galaxy( redshift=1.0, bulge=al.lp.Exponential( centre=(0.00, 0.00), ell_comps=(0.05, 0.05), intensity=1.2, effective_radius=0.1, ), mass=al.mp.Isothermal( centre=(0.0, 0.0), ell_comps=(0.05, 0.05), einstein_radius=0.6 ), ) source_galaxy = al.Galaxy( redshift=2.0, bulge=al.lp.Sersic( centre=(0.0, 0.0), ell_comps=(0.0, 0.111111), intensity=0.7, effective_radius=0.1, sersic_index=1.5, ), ) tracer = al.Tracer.from_galaxies(galaxies=[lens_galaxy_0, lens_galaxy_1, source_galaxy]) tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid_2d) tracer_plotter.figures_2d(image=True)
Chapter 1 of the HowToLens Jupyter notebook tutorials gives new users a full run through of lensing calculations with PyAutoLens and introduces those not familiar with gravitational lensing to the core theory.
A more detailed overview of these calculations can also be found in
Lens modeling is the process of taking data of a strong lens (e.g. imaging data from the Hubble Space Telescope or
interferometer data from ALMA) and fitting it with a lens model, to determine the
that best represent the observed strong lens.
Lens modeling with PyAutoLens uses the probabilistic programming language PyAutoFit, an open-source project that allows complex model fitting techniques to be straightforwardly integrated into scientific modeling software. Check it out if you are interested in developing your own software to perform advanced model-fitting!
We import PyAutoFit separately to PyAutoLens
import autofit as af
In this example, we consider Hubble Space Telescope imaging of the strong lens SLACS2303+1422. First, lets load this imaging dataset and plot it.
Note that the luminous emission of the foreground lens galaxy has been pre-subtracted in this image, to better highlight the source.
dataset_path = path.join("dataset", "slacs", "slacs2303+1422") imaging = al.Imaging.from_fits( image_path=path.join(dataset_path, "image.fits"), psf_path=path.join(dataset_path, "psf.fits"), noise_map_path=path.join(dataset_path, "noise_map.fits"), pixel_scales=0.05, ) imaging_plotter = aplt.ImagingPlotter(imaging=imaging) imaging_plotter.figures_2d(image=True)
We next mask the dataset, to remove the exterior regions of the image that do not contain emission from the lensed source galaxy.
Note how when we plot the
Imaging below, the figure now zooms into the masked region.
mask = al.Mask2D.circular( shape_native=imaging.shape_native, pixel_scales=imaging.pixel_scales, radius=3.0 ) imaging = imaging.apply_mask(mask=mask) imaging_plotter = aplt.ImagingPlotter(imaging=imaging) imaging_plotter.figures_2d(image=True)
2022-02-13 14:02:39,127 - autoarray.dataset.imaging - INFO - IMAGING - Data masked, contains a total of 11281 image-pixels
We now compose the lens model that we fit to the data using
Model objects. This behave analogously to the
MassProfile objects above, however their parameters are not specified and are instead determined
by a fitting procedure.
We will fit our strong lens data with two galaxies:
MassProfilerepresenting its mass.
LightProfilerepresenting a disk.
The redshifts of the lens (z=0.5) and source(z=1.0) are fixed.
lens_galaxy_model = af.Model( al.Galaxy, redshift=0.5, mass=al.mp.Isothermal ) source_galaxy_model = af.Model(al.Galaxy, redshift=1.0, disk=al.lp.Exponential)
We combine the lens and source model galaxies above into a
Collection, which is the model we will fit. Note how
we could easily extend this object to compose highly complex models containing many galaxies.
model = af.Collection(galaxies=af.Collection(lens=lens_galaxy_model, source=source_galaxy_model))
By printing the
Model's we see that each parameters has a prior associated with it, which is used by the
model-fitting procedure to fit the model.
print(lens_galaxy_model) print() print(source_galaxy_model)
Galaxy (centre_0, GaussianPrior, mean = 0.0, sigma = 0.1), (centre_1, GaussianPrior, mean = 0.0, sigma = 0.1), (ell_comps_0, GaussianPrior, mean = 0.0, sigma = 0.3), (ell_comps_1, GaussianPrior, mean = 0.0, sigma = 0.3), (einstein_radius, UniformPrior, lower_limit = 0.0, upper_limit = 8.0) Galaxy (centre_0, GaussianPrior, mean = 0.0, sigma = 0.3), (centre_1, GaussianPrior, mean = 0.0, sigma = 0.3), (ell_comps_0, GaussianPrior, mean = 0.0, sigma = 0.5), (ell_comps_1, GaussianPrior, mean = 0.0, sigma = 0.5), (intensity, LogUniformPrior, lower_limit = 1e-06, upper_limit = 1000000.0), (effective_radius, UniformPrior, lower_limit = 0.0, upper_limit = 30.0)
info attribute shows the model information in a more readable format:
We now choose the 'non-linear search', which is the fitting method used to determine the set of
MassProfile parameters that best-fit our data.
In this example we use dynesty, a nested sampling algorithm that in our experience has proven very effective at lens modeling.
search = af.DynestyStatic(name="overview_example")
2022-02-13 14:02:39,604 - autofit.non_linear.abstract_search - INFO - Creating search
To perform the model-fit, we create an
AnalysisImaging object which contains the
log likelihood function that the
non-linear search calls to fit the lens model to the data.
analysis = al.AnalysisImaging(dataset=imaging)
To perform the model-fit we pass the model and analysis to the search's fit method. This will output results (e.g., dynesty samples, model parameters, visualization) to hard-disk.
However, the lens modeling of this system takes 10-20 minutes. Therefore, to save time in this introduction, we have
commented out the
fit function below so you can skip through to the next section of the notebook. Feel free to
uncomment the code and run the lens modeling yourself!
Once a model-fit is running, PyAutoLens outputs the results of the search to hard-disk on-the-fly. This includes lens model parameter estimates with errors non-linear samples and the visualization of the best-fit lens model inferred by the search so far.
# result = search.fit(model=model, analysis=analysis)
The animation below shows a slide-show of the lens modeling procedure. Many lens models are fitted to the data over and over, gradually improving the quality of the fit to the data and looking more and more like the observed image.
Credit: Amy Etherington
The fit returns a
Result object, which contains the best-fit
FitImaging as well as the full
posterior information of the non-linear search, including all parameter samples, log likelihood values and tools to
compute the errors on the lens model.
Chapter 2 of the HowToLens Jupyter notebook tutorials gives new users a full run through of lens modeling, including a primer on Bayesian statistics, how a non-linear search works and effective strategies for fitting data of a strong lens.
More examples can be found in the
PyAutoLens provides tool for simulating strong lens data-sets, which can be used to test lens modeling pipelines and train neural networks to recognise and analyse images of strong lenses.
Simulating strong lens images uses a
SimulatorImaging object, which models the process that an instrument like the
Hubble Space Telescope goes through observe a strong lens. This includes accounting for the exposure time to
determine the signal-to-noise of the data, accounting for diffraction of the observed light by the telescope optics
and accounting for the background sky in the exposure which acts as a source of noise.
psf_2d = al.Kernel2D.from_gaussian(shape_native=(11, 11), sigma=0.1, pixel_scales=0.05) simulator = al.SimulatorImaging( exposure_time=300.0, background_sky_level=1.0, psf=psf_2d, add_poisson_noise=True )
Once we have a simulator, we can use it to create an imaging dataset which consists of an image, noise-map and Point Spread Function (PSF) by passing it a tracer and grid.
This uses the tracer above to create the image of the strong lens and then add the effects that occur during data acquisition.
imaging = simulator.via_tracer_from(tracer=tracer, grid=grid_2d)
By plotting a subplot of the
Imaging dataset, we can see this object includes the observed image of the strong lens
(which has had noise and other instrumental effects added to it) as well as a noise-map and PSF:
imaging_plotter = aplt.ImagingPlotter(imaging=imaging) imaging_plotter.subplot_imaging()