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 matplotlib.pyplot as plt
from os import path
import autolens as al
import autolens.plot as aplt
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 = al.Grid2D.uniform(
shape_native=(150, 150),
pixel_scales=0.05, # <- The pixel-scale describes the conversion from pixel units to arc-seconds.
)
grid_plotter = aplt.Grid2DPlotter(grid=grid)
grid_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
uses 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 Grid2D
and
create an image of the LightProfile
.
image = sersic_light_profile.image_2d_from(grid=grid)
plt.imshow(image.native) # The use of 'native' is described at the start of the HowToLens tutorials.
<matplotlib.image.AxesImage at 0x242fa364a30>
The PyAutoLens plot module provides methods for plotting objects and their properties, like the image of
a LightProfile
.
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
)
light_profile_plotter.figures_2d(image=True)
PyAutoLens uses MassProfile
objects to represent a galaxy's mass distribution and perform ray-tracing
calculations.
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)
The deflection angles are easily plotted using the PyAutoLens plot module.
mass_profile_plotter = aplt.MassProfilePlotter(
mass_profile=isothermal_mass_profile, grid=grid
)
mass_profile_plotter.figures_2d(
deflections_y=True, deflections_x=True
)
Many other lensing quantities are easily plotted with the MassProfilePltoter
, for example the convergence
and magnification
.
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 LightProfile
and 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)
lens_galaxy_plotter.figures_2d(image=True, deflections_y=True, deflections_x=True)
source_galaxy_plotter = aplt.GalaxyPlotter(galaxy=source_galaxy, grid=grid)
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 Galaxy
objects.
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 = tracer.image_2d_from(grid=grid)
tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid)
tracer_plotter.figures_2d(image=True)
The 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 Profile
's. A Tracer
'
can take many Galaxy
's.
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 lmp
module.
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)
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 notebooks/overview/overview_1_lensing.ipynb
.
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 LightProfile
's and MassProfile
's
that best represent the observed strong lens.
Lens modeling 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")
dataset = al.Imaging.from_fits(
data_path=path.join(dataset_path, "data.fits"),
psf_path=path.join(dataset_path, "psf.fits"),
noise_map_path=path.join(dataset_path, "noise_map.fits"),
pixel_scales=0.05,
)
dataset_plotter = aplt.ImagingPlotter(dataset=dataset)
dataset_plotter.figures_2d(data=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=dataset.shape_native,
pixel_scales=dataset.pixel_scales,
radius=3.0
)
dataset = dataset.apply_mask(mask=mask)
dataset_plotter = aplt.ImagingPlotter(dataset=dataset)
dataset_plotter.figures_2d(data=True)
2023-06-14 13:50:18,197 - autoarray.dataset.imaging.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. These behave analogously to the Galaxy
,
LightProfile
and 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:
Sersic
LightProfile
representing its light and an Isothermal
MassProfile
representing its mass.Exponential
LightProfile
representing a disk.The redshifts of the lens (z=0.155) and source(z=0.517) are fixed.
lens = af.Model(
al.Galaxy,
redshift=0.155,
bulge=al.lp.Sersic,
mass=al.mp.Isothermal
)
source = af.Model(al.Galaxy, redshift=0.517, 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, source=source))
By printing the Model
's we see that each parameter has a prior associated with it, which is used by the
model-fitting procedure to fit the model.
print(lens)
print()
print(source)
Galaxy (centre_0, GaussianPrior [12], mean = 0.0, sigma = 0.3), (centre_1, GaussianPrior [13], mean = 0.0, sigma = 0.3), (ell_comps_0, GaussianPrior [14], mean = 0.0, sigma = 0.3), (ell_comps_1, GaussianPrior [15], mean = 0.0, sigma = 0.3), (intensity, LogUniformPrior [16], lower_limit = 1e-06, upper_limit = 1000000.0), (effective_radius, UniformPrior [17], lower_limit = 0.0, upper_limit = 30.0), (sersic_index, UniformPrior [18], lower_limit = 0.8, upper_limit = 5.0), (centre_0, GaussianPrior [19], mean = 0.0, sigma = 0.1), (centre_1, GaussianPrior [20], mean = 0.0, sigma = 0.1), (ell_comps_0, GaussianPrior [21], mean = 0.0, sigma = 0.3), (ell_comps_1, GaussianPrior [22], mean = 0.0, sigma = 0.3), (einstein_radius, UniformPrior [23], lower_limit = 0.0, upper_limit = 8.0) Galaxy (centre_0, GaussianPrior [24], mean = 0.0, sigma = 0.3), (centre_1, GaussianPrior [25], mean = 0.0, sigma = 0.3), (ell_comps_0, GaussianPrior [26], mean = 0.0, sigma = 0.3), (ell_comps_1, GaussianPrior [27], mean = 0.0, sigma = 0.3), (intensity, LogUniformPrior [28], lower_limit = 1e-06, upper_limit = 1000000.0), (effective_radius, UniformPrior [29], lower_limit = 0.0, upper_limit = 30.0)
The info
attribute shows the model information in a more readable format:
print(model.info)
Total Free Parameters = 18 model Collection (N=18) galaxies Collection (N=18) lens Galaxy (N=12) bulge Sersic (N=7) mass Isothermal (N=5) source Galaxy (N=6) disk Exponential (N=6) galaxies lens redshift 0.155 bulge centre centre_0 GaussianPrior [12], mean = 0.0, sigma = 0.3 centre_1 GaussianPrior [13], mean = 0.0, sigma = 0.3 ell_comps ell_comps_0 GaussianPrior [14], mean = 0.0, sigma = 0.3 ell_comps_1 GaussianPrior [15], mean = 0.0, sigma = 0.3 intensity LogUniformPrior [16], lower_limit = 1e-06, upper_limit = 1000000.0 effective_radius UniformPrior [17], lower_limit = 0.0, upper_limit = 30.0 sersic_index UniformPrior [18], lower_limit = 0.8, upper_limit = 5.0 mass centre centre_0 GaussianPrior [19], mean = 0.0, sigma = 0.1 centre_1 GaussianPrior [20], mean = 0.0, sigma = 0.1 ell_comps ell_comps_0 GaussianPrior [21], mean = 0.0, sigma = 0.3 ell_comps_1 GaussianPrior [22], mean = 0.0, sigma = 0.3 einstein_radius UniformPrior [23], lower_limit = 0.0, upper_limit = 8.0 source redshift 0.517 disk centre centre_0 GaussianPrior [24], mean = 0.0, sigma = 0.3 centre_1 GaussianPrior [25], mean = 0.0, sigma = 0.3 ell_comps ell_comps_0 GaussianPrior [26], mean = 0.0, sigma = 0.3 ell_comps_1 GaussianPrior [27], mean = 0.0, sigma = 0.3 intensity LogUniformPrior [28], lower_limit = 1e-06, upper_limit = 1000000.0 effective_radius UniformPrior [29], lower_limit = 0.0, upper_limit = 30.0
We now choose the 'non-linear search', which is the fitting method used to determine the set of LightProfile
and MassProfile
parameters that best-fit our data.
In this example we use nautilus, a nested sampling algorithm that in our experience has proven very effective at lens modeling.
search = af.Nautilus(name="overview_example")
2023-06-14 13:53:29,570 - 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=dataset)
To perform the model-fit we pass the model and analysis to the search's fit method. This will output results (e.g., Nautilus samples, model parameters, visualization) to hard-disk.
However, the lens modeling of this system takes over 30 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 for this strong lens. The lens light has been subtracted beforehand, so that the animation clearly showed the fitting of the lensed source galaxy (e.g. the arcs).
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 Tracer
and 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 notebooks/overview/modeling.py
and notebooks/imaging/modeling
packages.
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 = 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, 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.
dataset = simulator.via_tracer_from(tracer=tracer, grid=grid)
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:
dataset_plotter = aplt.ImagingPlotter(dataset=dataset)
dataset_plotter.subplot_dataset()
The folder notebooks/imaging/simulators
includes many example notebooks for simulating strong lenses with a range
of different physical properties as well as creating imaging datasets representative of a variety of telescopes
(e.g. Hubble, Euclid).
Below, we show what a strong lens observation looks like for the lowest resolution instrument in these examples (the Vera Rubin Observatory) and highest resolution instrument (Keck Adaptive Optics).
Pixelizations reconstruct the source galaxy's light on a pixel-grid. Unlike LightProfile
's, they are able to
reconstruct the light of non-symmetric, irregular and clumpy sources.
The image below shows a pixelized source reconstruction of the strong lens SLACS1430+4105, where the source is reconstructed on a Voronoi mesh adapted to the source morphology, revealing it to be a grand-design face on spiral galaxy:
A complete overview of pixelized source reconstructions can be found
at notebooks/overview/overview_5_pixelizations.ipynb
.
Chapter 4 of the HowToLens lectures describes pixelizations in detail and teaches users how they can be used to perform lens modeling.
PyAutoLens supports modeling of interferometer data from submillimeter and radio observatories like ALMA and LOFAR.
Visibilities data is fitted directly in the uv-plane, circumventing issues that arise when fitting a dirty image such as correlated noise. This uses the non-uniform fast fourier transform algorithm PyNUFFT to efficiently map the lens model images to the uv-plane.
Given the irregular and clumpy nature of submm / radio sources, pixelized source reconstructions are vital for modeling interferometer datasets. This would be a slow process, however PyAutoLens* uses the linear algebra library PyLops to ensure visibilities fitting is efficient, even for datasets consisting of tens of millions of visibilities.
An overview of interferometer analysis is given in notebooks/overview/overview_6_interferometer.ipynb
and
the autolens_workspace/*/interferometer
package has example scripts for simulating datasets and lens modeling.
PyAutoLens supports multi-wavelength modeling of imaging datasets observed at different colors and combining these with interferometer data.
The appearance of the lens and source in a strong lens changes as a function of wavelength, therefore multi-wavelength lens modeling makes it easier to deblend their emission and offers more information to constrain the lens mass model.
An overview of multi-wavelength analysis is given in notebooks/overview/overview_7_mutli_wavelength.ipynb
and
the autolens_workspace/*/multi
package has example scripts for simulating datasets and lens modeling.
There are many lenses where the background source is not extended but is instead a point-source, for example strongly lensed quasars and supernovae.
For these objects, we do not want to model the source using a light profile, which implicitly assumes an extended surface brightness distribution. Instead, we assume that our source is a point source with a centre (y,x).
An overview of point-source analysis is given in notebooks/overview/overview_8_point_sources.ipynb
and
the autolens_workspace/*/point_source
package has example scripts for simulating datasets and lens modeling.
The strong lenses we've discussed so far have just a single lens galaxy responsible for the lensing. Group-scale strong lenses are systems where there two or more lens galaxies deflecting one or more background sources:
PyAutoLens has built in tools for modeling group-scale lenses, with no limit on the number of lens and source galaxies!
Overviews of group and analysis are given in notebooks/overview/overview_9_groups.ipynb
The autolens_workspace/*/group
package has example scripts for simulating datasets and lens modeling.
Congratulations, you've completed the PyAutoLens overview!
We recommend that new users first read through the overview examples of interest (autolens_workspace/notebooks/overview
)
and then find a start_here.ipynb
notebook for the feature they are interested in (for example, if you are looking to model
CCD imaging data, go to autolens_workspace/notebooks/imaging/modeling/start_here.ipynb
).
For users less experience with gravitational lensing, model fitting, Bayesian inference and scientific analysis you may wish to read through the HowToLens lectures. These teach you the basic principles of gravitational lensing and model fitting, with the target audience undergraduate students, masters students and above.
You can install PyAutoLens on your system and clone the autolens_workspace
tutorials
following the instructions on our readthedocs:
https://pyautolens.readthedocs.io/en/latest/installation/overview.html
Alternatively, you can begin using PyAutoLens in a web browser at the following Binder
link https://mybinder.org/v2/gh/Jammy2211/autolens_workspace/HEAD
.
A full overview of the autolens_workspace
folder structure is as follows:
notebooks/overview
: A more detailed overview of PyAutoLens's features .notebooks/howtolens
: A series of lectures aimed at undergraduates and above that teach gravitational lensing and
model fitting.notebooks/imaging/modeling
: Example scripts for fitting a lens model to CCD imaging data (e.g. HST).notebooks/imaging/data_preparation
A Preprocessing guide for preparing your CCD dataset for PyAutoLens.notebooks/imaging/simulators
: Simulating CCD imaging data.notebooks/imaging/results
: Tutorials explaining how to use the Result
object returned after lens modeling.notebooks/interferometer
: Interferometer modeling and simulations.notebooks/point_source
: Point-source modeling and simulations.notebooks/group
: Group-scale lens modeling and simulations.notebooks/plot
: Plotting examples and API reference .Examples describing advanced PyAutoLens features are also located on the autolens_workspace
(many are
in packages named advanced
). We advise that new users omit these packages until familiar with the software:
notebooks/multi
: Multi-wavelength modeling and simulations.notebooks/*/advanced/database
: Database tools for loading and analysing the results of large-scale lens model fits.notebooks/*/advanced/chaining
: Advanced modeling scripts which chain together multiple non-linear searches, including the SLaM (Source, Light and Mass) lens modeling pipelines.notebooks/*/advanced/log_likelihood_function
: A step-by-step guide to the log likelihood functions used by PyAutoLens.notebooks/*/advanced/graphcial
: Fitting graphical models to multiple imaging strong lens datasets simultaneously.notebooks/*/advanced/subhalo
: Dark matter subhalo modeling and simulations.